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DESIGN OF COMPUTATIONAL ALGORITHMS FOR 
OPTIMAL CONTROL BY HILBERT SPACE METHODS 

SUMMARY 


A new classification of computational methods for optimal control 
problems is given that includes the first variation methods of steepest descent, 
quasilinearization, boundary condition iteration, and the second variation 
method. The classification is based on viewing the solutions as a sequence 
of approximate control functions in the control space U, a sequence of approxi- 
mate adjoint functions corresponding to linear functionals in the dual space 
X* , or a sequence of approximate state functions in the state space, X. 

Three new computational algorithms for optimal control problems are 
developed, each belonging to one of the above categories. The first is based 
on Approximation in X* and is an extension of the quasi-Newton methods to 
boundary condition iteration. In contrast to former second-order methods 
for boundary condition iteration involving the Newton-Raphson method, the 
inverse Jacobian matrix of second partial derivatives is obtained by iteration. 
Only first partial derivatives of the cost functional are required, and a 
detailed technique far obtaining these derivatives is given. 

The second is based on Approximation in U. Concepts of asymptotic 
stability and Hamilton-Jacobi theory are used to obtain a sequence of func- 
tionals convergent to the optimal cost and which specify a minimizing sequence 
of control functions for the cost functional. A detailed algorithm for the state 
regulator problem is given. The control sequence is obtained by solving linear 
equations instead of nonlinear Riccati equations and possesses both monotonic 
and quadratic convergence. Numerical studies indicate that the algorithm is 
preferable to existing techniques, including Runge-Kutta integration of the 
Riccati equations and the Automatic Synthesis Program (ASP). 

The third is based on Approximation in X and exploits the concept of 
an inverse mapping of X into U. A second variation algorithm is formulated, 
and a detailed solution for the resulting boundary value problem is given. 

As an important special case the state regulator problem is treated, and an 
approximation to the optimal control is obtained. Comparison with the existing 
techniques involving the second variation indicates reduced computational 
requirements for problems in which the Euclidean dimension of U is less than 
that of X. 



A number of computational examples demonstrate the application of 
the new algorithms, and comparisons are made with the exact solutions. 

Finally, we obtain the solution to a class of state regulator problems. 
These problems involve first-order time-varying systems and are designed 
for testing computational algorithms. 
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CHAPTER 1 


INTRODUCTION 

Statement of Problem and Research Objectives 

The development of the Maximum Principle by Pontryagin and his 
colleagues in the 1950 r s is acknowledged as one of the most significant results 
in modern control theory. The method is a culmination of research into the 
Calculus of Variations by Bliss and the n Chicago school" more than 20 years 
earlier, and provides necessary (and often sufficient) conditions for an 
optimum. During the past 10 years, considerable effort has been devoted to 
exploiting the Maximum Principle for the design of computational algorithms 
for the numerical solution of control problems. 

A more fundamental approach is through the application of "functional 
analysis" and the concept of minimizing a linear functional on a normed linear 
vector space. Pontryagin, in fact, used these tools in the proof of his cele- 
brated result. The abstraction of the vector space approach yields greater 
insight into the problem structure, and the geometric character of the results 
provides a unified framework for further extensions. 

Basically the optimal control problem consists of the following items: 
a cost functional J to be minimized; dynamical system constraints described 
by differential, difference, or integral equations relating the state x and contol 
u; and constraints on the domain and range of these functions. A control u is 
said to be optimal if J is minimized and the constraints are satisfied. Through- 
out this report we assume that x and u belong to bounded open sets X and U, 
respectively, which are subsets of infinite dimensional Hilbert space. The 
system dynamics will be modeled by nonlinear ordinary differential equations. 
Terminal equality constraints on the state and unspecified terminal time will 
be considered. 

A computational solution to the optimal control problem stated above 
may be obtained by the following alternate viewpoints: (1) Approximation in 

U — choosing a sequence of approximate control functions in U by expressing 
the state in X as a bounded linear mapping of U into X and employing adjoint 
functions (representing elements inX*) to determine the mapping, (2) Approxi- 
mation in X* — choosing a sequence of approximate adjoint functions which 
correspond to linear functionals in X* , (3) Approximation in X — choosing a 
sequence of approximate state functions in X, provided that a mapping of X 
into U can be found. 
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The first objective of this research is to classify previous computational 
methods for optimal control problems according to the above scheme. The 
second objective is to formulate the optimal control problem as Approximation 
in X by carefi development of an "inverse" mapping of X into U. The last 
objective of the study is to develop new algorithms in each of the three cate- 
gories and demonstrate their application by computational examples. 

Organization of the Report 

In Chapter 1, the problem and research objectives are stated, and 
the organization of the report is presented. Basic terminology and definitions 
are introduced. 

In Chapter 2, the optimal control problem is defined, and assumptions 
are given. Using basic concepts from the theory of linear vector spaces 
(linear operators, Frechet differentials, adjoints, etc.) necessary conditions 
for a minimum of J on X x U x T are derived. As an example, the state 
regulator problem is reviewed, and known results are obtained and discussed. 

In Chapter 3, the computational solution of the optimal control problem 
is treated as Approximation in X* and U. These two categories are treated 
together because of their natural dependence on the necessary conditions for 
an optimal control. Well-known methods of descent for finite dimensional 
minimization are reviewed as a prototype for methods in infinite dimensions 
which follow, including the methods of steepest descent, boundary condition 
iteration, quasilinearization, and second variation. A computational algorithm 
for Approximation in X* is formulated based on boundary condition iteration by 
a quasi- Newton search. Next we describe a method for Approximation in U 
based on concepts of asymptotic stability and Hamilton-Jacobi theory. As an 
important special case of the latter method, a computational algorithm for the 
state regulator problem is developed. 

In Chapter 4, the inverse mapping from X into U is defined, and 
several important properties of the mapping are proved. A second variation 
method is developed based on the inverse mapping. The basic second 
variation algorithm is formulated and, as a special case, the state regulator 
problem is treated. 

In Chapters 3 and 4, the application of each new algorithm is demon- 
strated by computational examples, and comparisons are made with the exact 
solutions. 

In Chapter 5, conclusions and recommendations for future research 
related to the report topic are given. 
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In Appendix A, methods are given for solving the linear differential 
equation associated with the algorithm for monotone Approximation in U . 

In Appendix B, methods are given for solving the Accessory Problem 
associated with the second variation algorithm for Approximation in X. 

In Appendix C, we present a class of state regulator problems with 
nontrivial solutions that may be used to check computational algorithms. 


Terminology 


For the purpose of introducing relevant terminology and notation, the 
following definitions from the theory of linear vector spaces [ 1 ] are provided: 

Euclidean space IR n is the space of real valued n-tuples 

/ “ \Va 

(x ls ...,xj, (yi,...,yj with norm || x H^l x - I induced by the 


n 


n 


n 

inner product lx, y I = Y x.y.. Hilbert space L 2 (T) is the (complete) 

^ • 1 i=l 

space of IR - valued functions x, y which are Lebesgue square integrable on 
the closed interval T = [ t 0 , t ± J of IR, the extended real line, with norm 


x 


/h n \y 2 ti n 

=1 f Yj x ?(t) dt] induced by the inner product fx, yl = f Y 

\to i=l 1 / J t 0 i=l 


x ^(t)yi(t)dt . Note that IR n is a special case of Hilbert space. 

The interval [t 0 , t^] is denoted by T. The terms mapping and operator 
will be used interchangeably. A functional is a mapping of Hilbert space into 
IR. A linear mapping K of Hilbert space X into Y is bounded if there exists a 
constant a such that ||Kx|| ^ ce 1 1 x 1 1 for all xeX . The dual space X* of 
Hilbert space X consists of all bounded linear functionals on X. The adjoint 
K* of a bounded linear mapping K : X — Y where X and Y are Hilbert space is 

a mapping K* : Y — X defined by £y, Kx J = y,xj for xeX , yeY . A 

bounded linear operator K is self adjoint if K = K* . A self adjoint matrix 

operator K on IR n is positive d e finite (positive semidefinite) if |x, KxJ > 0(^0) 

for all x + , KxJ = 0 for x = 9 ^ and will be denoted by K > 0(^0) . 


.n 


A self adjoint matrix operator on IR is called symmetric. The transpose of a 
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matrix operator M on IR n is denoted by M* in accordance with the definition 
of adjoint. The null operator is denoted 0 or 6 and identity by I where the 
space is obvious by the context. 


Differentiation of a function x on T with respect to the real variable t 
will be denoted by x. Partial derivatives will be denoted by subscripts as in 
the following example: 


L 

x 


8L 


9L 


9x t 9x 


,L 


n 


xx 


9 2 l 

axjSxj 

9 2 L 

9x 9Xi 
n 1 


9 2 L 

9xj9x n 

9 2 L 

9x 9x 

n n 


where L : IR n — IR. The inverse of matrix M is denoted by M -1 . If X and 
Y are subsets of Hilbert space, the cartesian product X x Y is the collection 
of ordered pairs (x,y) with xeX , yeY . Addition and scalar multiplication 
are defined by (x^yj) + (x 2 ,y 2 ) = (x* + x 2 ,y 1 + y 2 ) and a(x,y) = (ax, ay). 
( k) 

C (T) is the class of continuously k-differen tiable functions on the interval 
T. Given X and Y Hilbert spaces, the first Gateaux differential (the first 
variation) 6F(x;h) of F : X — Y at x with increment h is defined by 


1 J 

<5F(x ; h) = lim — (F(x + ah) - F(x) ) = — F(x + cxh) 
„ a da 

a-* 0 


a = 0 


if the limit exists in the sense of convergence in the norm. If for fixed xeX 
and each heX there exists <5F(x;h) which is linear and bounded and 


lim 


II h 


(l|F(x +h) - F(x) - <5F(x, h) ||) = 0 


then F is said to be first Frechet differentiable at x, and the unique mapping 
dF(x^h) = 6F(x;h) is called the first Frechet differential of F at x with 
increment h. Higher-order differentials are defined in a similar manner. We 
shall say that F has a relative weak (strong) minimum at xeX if there exists 
an a > 0 such that F(x) > F(x) whenever |Jx - x|| < a and xGC 1 (T) 

(x e C° (T) ) . For brevity we shall often say that F has a minimum over 
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all X if the above condition holds for all xGX . A sequence of functions 

( ki 

x , k = 0,1,... . in Hilbert space X is called a minimizing sequence 
for the functional V : X -*• IR if 


lim v(x (k) ) 

k — oo 


= inf V (x) 
x€X 
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CHAPTER 2 


ANALYSIS OF THE OPTIMAL CONTROL 
PROBLEM IN HILBERT SPACE 


In this chapter, the optimal control problem is defined. Using elemen- 
tary concepts from the theory of linear vector spaces, necessary conditions 
for a minimum of J on X x U x T are derived. As an example, the state 
regulator problem is reviewed, and known results are obtained and analyzed. 


Definition of the Optimal Control Problem 

Consider the problem of determining a control function u in the 
control space U = L^(T), T = [t 0 , tj such that a cost functional 

J:XxUxT — IR defined by 


J(u,x) = f L(x,u,t) dt 


( 1 ) 


to 

is minimized for all x in the state space X = L 2 (T) , subject to nonlinear 
dynamical system constraints 

x = f (x , u , t) , x (t 0 ) = c (2) 

where f:X x U x T — X and subject to nonlinear terminal equality constraints 

ip(x , t l ) = d , (3) 

where il : X x T — IR^. Restrictions on the problem are considered in the 
following section. 


Basic Assumptions 

The basic assumptions are: 

1. The pair of admissible functions (u, x) and their increments 
(6u, 6x) that satisfy the linearized dynamical system 
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( 4 ) 


8x = f (x , u , t) Sx + f (x , u , t) <5u , 6x(t 0 ) = 9 

X u 

are continuously differentiable on T . 

2. Mappings L, f, and ip are continuously second differentiable on 
their domains. 

3. The functional J possesses first and second Frechet differentials 
and is a convex function of its arguments over all X x U x T. 

4. The linearized dynamical system [eq. (4)] is uniformly completely 
controllable on T [2]. This condition is equivalent to requiring that the 
controllability matrix 

h 

f $ (to , s) f (x , u , s) f * (x , u , s) 4>' ,s (t 0 , s) ds (5) 

J u u 

is positive definite where <f> (t, s) is the transition matrix of equation (4) and 
satisfies the matrix differential equation 

4r $ (t , t 0 ) = f (x , u , t) $ (t , t 0 ) , *<to,to)=I . (6) 

ut x 

5. The terminal equality constraints are linearly independent. Thus 

rank |^ x (x,tj)| =p . (7) 

Condition 3 is a sufficient condition for a minimum of J and may be 
too restrictive for applications. Furthermore, conditions 1 and 2 may also 
be too strong and should be examined in specific cases. 

If conditions 1 through 5 are satisfied, the following development 
provides both necessary and sufficient conditions for a weak minimum of J on 
X x U x T . 


Necessary Conditions 

To facilitate the presentation of the concepts to follow, let us assume 
that the terminal time t 1 is specified in advance. In general, it is only 


9 


available implicitly through the terminal equality constraint [eq. (3)]. 
Necessary conditions for the more general case have been given in References 
3 and 4. 

Let the Lagrangian functional F : X x U x X* x IR P x T — IR be 
defined by 

F (x , u , A , v , t) = J (x , u) + £a , f (x , u , t) - x J 

+ [p , ip (x, t t )J^ , (8) 

where xeX = L^(T) and ueU = L^(T) . Elements, A, v are the unique 
"Lagrange multiplier" or adjoint functions, which represent elements in the 
dual spaces X* and (IR P )* , respectively. 1 

Suppose we define the bounded mapping, K : X -*X, for all xeX 
by 


t 

(Kz) (t) = J z(s) ds . (9) 

to 

Then the state function may be expressed in terms of its derivative as 

x(t) = c + (Kx) (t)l (10) 

where c is the initial condition for the dynamical system [eq. (2)J. Equation 
(10) implies that F as defined by equation (8) may be viewed as a function of 

x. 


If equations (2) and (3) are satisfied, a necessary condition for a 
weak minimum of F , and consequently a minimum of J over the product 
space X x U x T , is that the first Frechet differential of F at x and u vanish. 


1. Because of the (isomorphic) equivalence of A, v with elements of 
their dual spaces [1], we shall write AeX* , v e (IR P )* and we shall 
identify X* and (IR P )* with L^(T) and IR P , respectively. 
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Existence of the Frechet differential implies the existence of the Gateaux 
differential [5] , and the latter provides a convenient means for computing 
dF . Hence, for all SxeX , 


dF(x ; <5x) = ~ f(c + K (x + adx) , u, X , v , t) j a ~ q 



‘ ir( [L x (x,u,t) , 6x] 

+ £ X , f (x , u , t) 6x ^ 

- Jx , x J-fx, a6x] +0(|| 6x || 2 ) 

+£ v , (x , tj) 6x(t I )J^ +0(|| 6x(ti) || 2 )^ j a = 0 


2 

where L 0 , f 0 , % are constants, and 0(|| 6x || 2 ) , 0(|| 6x(t 1 )|| 1 ) are 
remainder terms for first-order Taylor series expansions of L , f , ip at x . 
By substituting KSx for 6x and differentiating with respect to a , we 
obtain for all 6x X , 
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dF(x ; Sx) = £ L x (x, u, t) , KSx J + £ A , iMx , u , t) K<5x J 
- Ja, <5x J + , ip^ (x, %) (Kdx)(tt)J 

= £ K* L x (x , u , t) + f* (x , u , t) A -A, 6x] 

+ J (x , tj) v , (K6x) (tj)J ^ 

= £ K* (x , u , t) + f^ : ' (x , u , t) A ) - A + (x , tj ) v <5x J 


= 9 


In computing dF(x, 6x) we have performed the following steps: (1) expanded 

L, f, ip in a first-order Taylor series about x; (2) dropped the constant terms 
L 0 , fo> which are independent of a; (3) used the Hilbert space inner pro- 

Jl 

duct notation [,] for J [,] dt ; (4) used equation (10) to write 6x as a func- 

to 

tion of its derivative 6x; (5) performed the indicated differentiation with respect 
to a ; (6) used the identity 

£ ^ (x , tj) v , (Kdx) (tj ) J j = J (x,tj)v , 6 x J . (13) 

Similarly, for all 6ueU , 
d 

dF(u;6u) = — F (x, u + adu, \ , t) = Q 

= [l u (x, u, t) , <5u J+J A, f u (x, u, t) 6uJ (14) 


= e 

A well-known theorem in the Calculus of Variations is the Euler-Lagrange 

Lemma, which states if /3 is a continuous function on T and if J/3, 6zJ = 0 

for every 6z in C 1 (T) where 6z(t 0 ) = 6z(t 1 ) = 9 , then /3 = 9 on T 
[6, 7] . This theorem is of fundamental importance and will be used frequently 
throughout this report. 
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Therefore, by the Euler- Lagrange Lemma, equations (12) and (14) 
imply that an optimal control function u satisfies 

L (x , u, t) + f * (x, u, t) A = e , (15) 

u u 

and the adjoint function A satisfies 

A = K* ( L (x , u , t) + f * (x , u , t) A ) + #* (x , % ) v . ( 16) 

The adjoint operator K* exists since K is bounded on X. The mapping 
K has been used by Mitter [8] to obtain necessary conditions and to develop 
an iterative procedure for determining the optimal control function. This 
section is based on Gruver [9] and is an extension of the results contained 
therein. A less formal treatment is given by Blum [10] . 


State Regulator Problem 

A well-known problem in optimal control concerns minimizing on 
X x U x T a quadratic cost functional 

J(u , x) = \ ^ £x , QxJ + ju , Rujj ^ (17) 

subject to linear dynamical system constraints 

x = Ax + Bu , x(t 0 ) = c , (18) 

where Q, R, A, B are real valued nxn, rxr, nxn, nxr matrix 
operators, respectively, on T . It is assumed that Q 0 and R > 0 and 
that the pair (A, B) is uniformly completely controllable on T . Hence, 
the conditions of the section, Basic Assumptions, are satisfied, and equations 
(15) and (16) are both necessary and sufficient for a unique minimum of J . 

From equation (15) the optimal control ft is 

u = -R -1 B?A . (19) 

Equations (16), (18), and (19) imply that x and A satisfy the integral 
equations 


t 

x (t) = c + f (Ax - BR _1 B* A) ds (20) 
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( 21 ) 


X(t) = / (Qx+A*X)ds 
t 

where elements of A, B, Q, x, X are functions of the real variable s. In 
the last equation the adjoint of the integral operator K in equation (9) has been 
computed [1, p. 154]. In 1960 Kalman showed that equations (20) and (21) 
may be decoupled by employing the "Riccati transformation" 

X = Px , (22) 

where P is an n x n real symmetric positive definit matrix operator on T. 

By substituting equations (20) and (21) into equation (22), differentiating 
equations (20) and (20) with respect to t, and eliminating X using equation 
(22) , the optimal control function may be expressed as 

u = -R _1 B*Px 

where P satisfies the matrix Riccati equation 

- P = PA + A* P - PBR _1 B* P +Q , P(t)=© . (23) 

Furthermore, the minimum cost functional has the explicit representation 
J (u , x) = i £x , Px] f 

and lim J(u, x) is a Lyapunov functional for the optimally controlled system 

tj — «• + OO 

( 2] . Suppose that A, B, Q, R are constant matrices on T and tj — + °° . 
Because of the time invariance, we may consider the problem on the infinte 
interval [0, ■»] by translation of the orgin. Then, Kalman has shown that 
"P = lim P(t) exists, is unique, and satisfies the matrix equation 

t — - =o 


PA + A* P - PBR -1 B* P + Q = © . (24) 

Furthermore, there exists a solution to equation (24) such that the controlled 
system is asymptotically stable. 

As an alternative to solving equation (23) , the adjoint function X can 
be obtained in terms of a 2n x 2n transition matrix f 1 1 ] . One difficulty is 
that the procedure requires inversion of submatrices of the transition matrix, 
which may become ill conditioned because of roundoff error in computation. 
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Kalman has also obtained an algorithm for computing discrete time optimal 
control functions by solving Riccati-type difference equations [12, 133. In 
the section, The State Regulator Problem Revisited, I, Chapter 3, we describe 
another technique that was developed by the author and that leads to a monotone 
and quadratic convergent approximation of the minimum cost functional and the 
optimal control for equations (17) and (18) by solving linear matrix differential 
equations instead of nonlinear Riccati equations. The method is applicable to 
both time-invariant and time-varying systems. 
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CHAPTER 3 


APPROXIMATION IN X’ AND U 


In this chapter, we consider the computational solution of optimal 
control problems first, as Approximation in X* and second, as Approxima- 
tion in U . These two categories are treated together because of their natural 
relationship to the necessary conditions [eqs. (15) and (16) J. Methods of 
descent for finite dimensional minimization are reviewed as prototypes for 
infinite dimensional spaces. As realizations of the latter, we examine the well- 
known "direct" methods of steepest descent and second variations and also 
the "indirect" methods of boundary condition iteration and quasilinearization. 

A second-order computational algorithm for Approximation in X' 1 ' 
by conjugate direction search is formulated. An accurate means for computing 
the gradient is given. Next, we describe a method for Approximation in U 
based on concepts of Hamilton- Jacobi theory and asympotic stability. As an 
important special case, a computational algorithm for the state regulator 
problem is developed that requires solution of linear matrix equations instead 
of nonlinear Riccati equations. Methods for solving these equations are 
discussed. 


Descent Algorithms for Euclidean Space 

In this section we offer a preview of the functional minimization prob- 
lem based on well-known methods for unconstrained minimization in finite 
dimensional Euclidean space. Suppose there is given a functional 

J : IR n — IR , and the problem is to find a minimum of J on IR n . It is 
assumed that the gradient g = J may be obtained analytically as a function. 

X 

Throughout this section the gradient, the IR valued function J defined by 

X 

dF (x: h) [A h ], , will be denoted by g . 

A basic consideration in the design of minimization algorithms is 

finding a systematic method for updating an element x^ i n lR n such that 

| x (k) J converges to a point where J is minimized. Essentially, the update 

consists of a step direction in IR n and a step length. In general, the step 
length must be determined by a one-dimensional search procedure such as 
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quadratic or cubic interpolation, Fibonacci search, or golden section minimi- 
zation [14, 15], Choice of the step direction is more complicated, and one 
classification scheme is as follows: (1) first-order methods; (2) second- 

order methods; (3) conjugate direction methods. In (1) the step direction 

in IR n for the update is restricted to the negative gradient of J (the direction 

of "steepest descent" in IR n ) and, hence, only first-order partial derivatives 
of J are needed. In (2) the step direction employs second-order partial 
derivatives of J . In (3) the step directions are chosen g - conjugate as 

X 

defined below and require only first-order partial derivatives of J . More- 
over, convergence of a minimizing sequence is quadratic in the sense of 

minimizing a quadratic functional on lR n in at most n iterations. 
FIRST-ORDER METHODS 

A model for first-order methods is steepest descent. In this technique 
the step direction is restricted to the negative gradient of J . The procedure 
may be implemented either in a continuous or a discrete manner. The latter 
is well-suited for digital computer calculation, and the steps in the algorithm 
are as follows: 

1. Choose an initial element, x^; let k = 0 . 

2. Compute the step direction, 


(k) , (k) . 

s = -g (x ) 


(25) 


(k) 


3. Compute (by a one-dimensional search) a step length a which 

(h 

a 

(k) 


'( 5 


... (k) (k) (k) 

minimizes J| x + a s 


■Ullllly 

) ' 


4. Update x 


according to 


(k+i) (k) (k) (k) 

x = x + a s 


(26) 


.(k), 


5. Let k — k + 1 and repeat steps 2 through 4 until ||g(x 7 ) | [ t is 

less than a predetermined positive number. 


The method of steepest descent is extremely stable in a large neighbor- 
hood of a minimum. However, computational results indicate slow convergence. 



except for the case in which the "level curves" (loci of constant J in IR n ) 
are hyperspheres. Since successive steps are orthogonal, the procedure often 
leads to inefficient "zig-zagging". Mathematically, it has been shown that 

convergence of j J i s at least as fast as a geometric series with ratio 

(M-m)/M where ml<g < MI , 0<m<M [16], 

X 

SECOND-ORDER METHODS 

A model for second-order methods is the classical Newton-Raphson 
method. Basically, it is an iterative technique for solving for the "roots" of 

an equation T (x) = 0 where T is a nonlinear mapping of IR n into IR n . 

Suppose J : IR n — IR , then the gradient equation, 


g(x) = 6 , (27) 

A 

is a necessary condition for a minimum of J at x , and the Newton-Raphson 
method is a convenient means for obtaining the minimizing element x. The 
algorithm is identical to that given for steepest descent except that steps 3 
and 4 are replaced by the following iteration: 


(k+i) 

x 




) g(x 



(28) 


where g is the Jacobian matrix of g . 2 Conditions for existence of the 

inverse Jacobian, and convergence of the minimizing sequence j x^ j are 

given in Reference 17. These conditions are rather lengthy, however, and in 
practice are usually not checked. The Newton-Raphson method is very efficient 

close to the minimum since g is a measure of curvature in IR n . In fact. 

x 

ii J is a quadratic functional on IR , then under certain conditions the 
algorithm converges in one iteration. A disadvantage, however, is that the 
inverse Jacobian may fail to exist during the descent, whereas if the problem 


2 . 


The n x n Jacobian matrix g of the IR n -valued function g = J is 

x x 

also called the Hessian of the quadratic functional J . 
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possesses a well-defined gradient, it could probably be solved by steepest 
descent. An additional disadvantage is that the iteration equation (28) may be 

unstable for certain choices of and fail to converge. This instability 

may sometimes be corrected by adjustment of the step length. 

CONJUGATE DIRECTION METHODS 


A compromise between the methods of steepest descent and Newton- 
Raphson is the following class of minimization techniques, which possesses 
quadratic convergence — in the sense of minimizing a quadratic functional 


n 


on IR in at most n steps. 


The class has the advantage that only first-order partial derivatives 
of J are needed. The earliest conjugate direction method, called the method 
of conjugate gradients, was developed by Hestenes and Stieffel [18], 

The name, "conjugate direction methods, " arises because successive 


directions s 


(k) 


k = 0,1,2, 




in IR are chosen to be conjugate 


(orthogonal with respect to the Jacobian matrix, g ) as defined by 


l 


(j) 

s . g l x 

A 


.0 


lk) ) s 


(k) 


] 


= 0 


j * k 


(29) 


As in the Newton-Raphson technique, an essential simplification occurs if J 


n 


is a quadratic functional on IR . In fact, if J is of the form 


J(x) = c + 


[»-],♦ l[ x ' Ax ], 


(30) 


_n 


where b e IR and A is a positive definite matrix, then the main theorem 

(k) 


for the method specifies a minimizing sequence { x 
property that / x^^ — A_1 


for J with the 

A -1 b [1, p. 294]. Convergence of the method is 

far superior to steepest descent. Details concerning the construction of the 
minimizing sequence, the step direction, and step length are available in the 
references and will not be needed in the sequel. 


Another important class of conjugate direction techniques is the quasi- 
Newton methods. Suppose J is a quadratic functional as in equation (30) . 
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Then this class of methods generates a sequence of positive definite n x n 

matrices j 1 } which converge to the inverse Jacobian A -1 . Many 

(k) 

different schemes for updating H have been proposed, although most 
computational results and the greatest success has been reported for that 
given by Davidon and later refined by Fletcher and Powell [19]. The authors 
of the latter reference have also proved the stability and convergence of the 
algorithm for nonlinear J . Details concerning the construction of the minimi- 
zing sequence are given in the references [19, 15], 

The method of conjugate gradients requires saving the gradient between 
iterations, and computer storage requirements are about the same as for the 
method of steepest descent. Quasi-Newton methods, however, provide faster 
convergence for nonlinear problems in return for storage between iterations 

(k) 

of both the gradient and the n x n matrix H . 

Descent Algorithms for Hilbert Space 

In the following section, computational methods for optimal control are 
classified according to (1) methods based on the first variation of the cost 
functional; (2) methods based on the second variation of the cost functional. 
Within these categories, particularly in (1), the first, second order, and 
conjugate direction methods of the section, Descent Algorithm for Euclidean 
Space, may be used to compute the control function update. 

METHODS BASED ON THE FIRST VARIATION 

We now consider minimizing a functional J:XxUxT — IR where 

a minimizing sequence j u^ j for J lies in a bounded open subset U of 

Hilbert space. Suppose we are given the optimal control problem of the section, 
Definition of the Optimal Control Problem, Chapter 2, of minimizing a cost 
functional 


J(u,x) = f L(x,u,t) dt 
to 

subject to nonlinear dynamical system constraints 
x = f(x , u , t) , x(t 0 ) = c 


(31) 


(32) 
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and terminal equality constraints 


ip(x , t i) = 0 . (33) 

In general the terminal time may be free and, consequently, determined 
implicitly through equation (33). 

Necessary conditions for a solution to a special case of this problem, 

(t t specified) were derived in the section, Necessary Conditions, Chapter 2. 
Suppose we define the Hamiltonian functional H : X x X v x U x T — IR by 

H(x , A , u , t) = L(x , u , t) + £ A , f(x , u , t)J (34) 

where AG X ' (see footnote 1 in Chapter 2) . Necessary conditions for a 
weak minimum of J over the product space X x U x T are derived in 
References 3 and 4, and the results are as follows: 


A 


x = H (x , A , u , t) 
A 

c+ 

o 

II 

O 

(35) 

• A 

A = -H ( x , A , u , t) 

x 


(36) 

CD 

II 

S3 

P^ 

>> 

p > 

c+- 


(37) 

A(tj) = ( x , tj) v 

X 


(38) 

A 

0 = H(x , A , u , tj) + ' 

^ t (x , tj) , 

(39) 

e = ipix , tj) 


(40) 


where v is a constant adjoint function in (IR^) " . Notice that if tj is 
specified, equations (39) and (40) are not needed, and equations (36) through 
(38) are equivalent to equations (15) and (16). Equations (35) through (40) 
constitute a two-point boundary problem since A (to) is unknown. Thus a 
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control function u e U is optimal provided that an adjoint function A can be 
found such that the following conditions are satisfied: 

1. The canonical equations (35) and (36). 

2. The min-H requirement, equation (37). 

3. The terminal boundary conditions, equations (38) through (40). 
Example 1 

As a concrete example of the necessary conditions just developed, con- 
sider the minimization of a quadratic cost functional 

1 jfi 

J(u , x) = — f (x 2 + u 2 ) dt 

2 0 


subject to a linear first-order dynamical constraint 


x = - x + -u , x( 0) =1 

and terminal equality constraint 

x( t j) -5 = 0 , 

where t t is free. For this simple problem, the Hamiltonian functional is 
H(x , X , u , t) = (x 2 + u 2 ) + X( -x + u) , 


and equations (35) through (40) are as follows: 
1. The canonical equations 

* 

x = -x + u , x ( 0) = 1 

A. = -x + A 


(35') 

(36') 
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2. The min-H requirement 


u + A = 0 . (37') 

3. The terminal boundary conditions 

A(t t ) = v (38') 

0 = ^-^-(x 2 + u 2 ) +A(-x + u) + xuj | ^ (39') 

0 = x(tj) - 5 . (40') 


In the past, an iterative solution to equations (35) through (40) has involved 
satisfying any two of the preceding conditions while adjusting the third. Thus we 
obtain the well-known methods of quasilinearization (adjust the adjoint functions 
in condition 1) £ 20 j , steepest descent (adjust the control function in condition 2) 
[21, 22] , and boundary condition iteration (adjust the adjoint functions in condi- 
tion 3) [23,24]. Quasilinearization and boundary condition iteration are called 

"indirect" methods since a minimizing sequence j u^ j for J is obtained 

indirectly by iteration of the canonical equations (35) and (36) , usually by 
means of a generalized Newton-Raphson search [17]. In this work we feel that 
Approximation in X v is a more natural classification since the basic iteration 
occurs in X* . On the other hand, steepest descent is referred to as a "direct" 
method since the Hamiltonian functional and, consequently, the cost functional 
are minimized directly in the control space U . Thus steepest descent may 
be classified as Approximation in U . 

Usually the adjustment in X* or U is enforced by a form of gradient 
descent. In some cases the canonical equations (35) and (36) can be expressed 
as a contraction mapping on the Hilbert space X x X" x U x T , and one 
could employ the method of successive approximations attributed to Picard and 
formalized by Banach to obtain a fixed point. Some results using this concept 
have been obtained by deJong [25]. 

Let us consider the boundary condition iteration method in more detail 
for the special case tj specified and unconstrained terminal state. The general 
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case is given in References 26 and 24. Suppose equation (37) may be solved 
explicitly for u as a function of x and A . Then, by substituting u(x, A) 
into equations (35) and (36), the canonical equations may also be expressed in 
terms of x and A as follows: 


H (x , A , u(x , A) , t) , 

A 

o 

II 

CD 

(41) 

-H (x , A , u(x , A) , t) , 

X 

Mt„) = d 

(42) 


where d is the initial condition for the adjoint function and is unknown. The 
basic boundary condition iteration algorithm for Approximation in X* pro- 
ceeds as follows: 

Choose an initial value for d^ ; let k = 0 . 

Integrate equations (41) and (42) forward from x(t 0 ) , A(t 0 ) . 

(k) (k) 

Compute the function u from equation (37) and store x v ' (t t ) , 

(k) 

Update d in the direction of steepest descent of the terminal 
boundary conditions. 

5. Let k — k + 1 , and repeat steps 2 through 4 until equations (38) 
through (40) are satisfied to within a predetermined bound. 

In step 4, we usually minimize a functional E which involves the terminal 
boundary conditions. A procedure for choosing this functional will be explained 
in Example 2 and, for a more general case, in the forthcoming section, A 
Conjugate Direct Method for Approximation in X" . 

Example 2 

Use of the boundary condition iteration algorithm will now be demon- 
strated by a special case of Example 1 given earlier in this section. Suppose 
that the terminal time is specified at tj = 1 , and the terminal state is uncon- 
strained. Then, equations (39') and (40') are not needed. From the min-H 
requirement, equation (37'), the control function may be expressed as u = -A 
and used to eliminate u in equations (35 1 ) and (36’). The result is 


1 . 

2 . 


4. 


x = -x - A , x( 0) = i 


( 41 ') 
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A 


-x + X 


A(0) = d 


(42') 


where d is the unknown initial condition for the adjoint function. A convenient 

1 , 

choice for a functional of terminal error is E(x, A,t) = —X (t) since if 

E(x,A,t 1 ) = 0, the terminal boundary condition is satisfied. To implement 
a descent algorithm, the gradient E^ = E + E^A^ is needed. Let us 

differentiate equations (41 1 ) and (42') with respect to d . The result is the 
"sensitivity equations," 



x ( 0) = 0 

d 


A d = 


= -x + A 
d d 


A ( 0) = 1 
d 


a set of coupled first-order differential equations that may be solved for A^ . 

A rigorous derivation of these equations has been given by Levine [23]. Suppose, 
for example, that d = 0 . Then by solving equations (41*) and (42') for the 
functions x , A and solving the "sensitivity equations" for jc . , X d , the 
gradient is 


E , = A A . 


-n/ 2 t \l~2 t 
e - e 


U-N/27e"^ t -(l+^e^ t 


Now that the preliminary equations have been set up, the first iteration of the 
boundary condition algorithm proceeds as follows: 


1. Let d = 0 . 


2. Integrate equations (41') and (42') forward from 
( 0 ) _ 1 


x 


2 42 


nrv ^2t ,, r- -\f2t 

( 1 - 42) e - ( 1 + 4 2) e 


( 0 ) 


2\f2 


42 t -42 t 
e - e 
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3. Compute u from equation (37'): 


u 


( 0 ) 


i / Vat -V2t 

e - e 


2 V2 


and store 


( 0 ) , .. 

x (1) 


— (-o. 

2\l~2 \ 


\[2 -42 

586 e - 2. 414 e 


2\T2 \ 


) 


4. Update d^ according to 


d U> = d W + c< W E d (x<°> , x (0) , 1 ) 

= 4 ° l ( e -^. e ^)( 0 . 586 e -^ + 2 . 144 e ^) 

where a ^ is a constant chosen to minimize E(x,A, 1). 

Step 4 involves a one-dimensional search for , which requires 

repeated evaluation of E(x,A, 1) using the functions and A ^ from 

equations (41*) and (42*) with initial conditions x^(0) = 1 , A^(0) = d^ 
Steps 2 through 4 are repeated until E(x^ , A^ , tj) = < a x 

where a x is a small positive number. Notice that integrating the canonical 
equations has given rise to terms of the form ^ , functions which increase 
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with t . This illustrates an inherent weakness of the boundary condition 
iteration method, because the technique is applicable only to those optimal 
control problems for which the adjoint equation (36) is stable when integrated 
forward on T , or at least the instabilities do not predominate. Fortunately, 
this includes a large class of trajectory optimization problems for which 
atmospheric effects are neglibible. The algorithm is particularly attractive 
since only storage of the terminal values of x(tj) and X (t ± > is required 
between iterations. In contrast, steepest descent methods require storing at 
least the entire adjoint function A on T . 

In the past, choice of the update in step 4 was based on minimizing a 
functional of terminal error (such as a quadratic) using steepest descent 
[23, 27] , conjugate gradient search [28] , or the Newton-Raphson method 
[26,24] to satisfy the terminal boundary conditions. In the section, A 
Conjugate Direction Method for Approximation in X' 1 ' , we will describe the 
application of quasi-Newton methods described in the section, Conjugate 
Direction Methods, to the boundary condition iteration technique. 

Finally, for comparison with the latter algorithm the basic steepest 
descent algorithm for Approximation in U is formulated as follows: 

1. Choose a nominal control function u^ ; let k = 0 . 

2. Integrate equation (35) forward from x(to) . Then compute A(t t ) , 

(k) 

and store the function x x ' . 

3. Integrate equation (36) backward from A(t 1 ) , and store the function 

A (k) . 

4. Update the control function in the direction of steepest descent of H 
according to 


(k+1) (k) (k) / (k) (k) (k) \ 

u = u -a H lx ,u > t ) 

where a ^ is chosen to minimize ||H^x^ , A^ , u^ k+ *^ , t^|| 2 . 

5. Let k — k + 1 and repeat steps 2 through 4 until 
/ (k) (k) (k) \ 

[ | H (x v , A v , u , t)|| 2 is less than a predetermined positive number. 
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Example 3 


Use of the steepest descent algorithm is now demonstrated by the 
problem introduced in Example 2. The first iteration proceeds as follows: 


1 . 


Let u 


( 0 ) 


0 . 


2. Integrate equation (35') forward from x(0) = 1: 


3. Integrate equation (36 1 ) backward 3 from X ( 1) = 0 . 

(0) _ l/t-i -(t+l)\ 

X 2 (e " 6 ) 

4. Update the control function according to 

(1) (0) (0) a (0) ( t-1 -(t+l)\ 

u = u - - a H = — - — ( e - e I , 

u 2 \ / 

where a ^ is chosen to minimize 

II H (x (0 ', X (0) , u ( °> , t)ll* - /(f ((x <0) ) 2 + (» (1) ) 2 ) 

+ X (0) (-x (W + u (1) )) ! d t 
using a one-dimensional search procedure [14, 15]. 


3. To integrate backward, we make the change of variable t = 1 - s in 

equation (36'). Therefore, dx/dt = - dX/ds , and equation (36') becomes 

, / , s-1 

dX/ds + X = e 
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Steps 2 through 4 of the algorithm are repeated until 

IIH^H 2 = J dt < cti 2 • It is significant that the adjoint 

equation is integrated in its "stable" direction. Consequently, the steepest 
descent method may be applied to a wider class of problems than the boundary 
condition iteration method in return for storage of at least the entire adjoint 
function X over the interval T . In terms of numerical integration using a 
digital computer, this requires discretizing the interval T and storing num- 
bers at many points. 

In step 4, the conjugate direction methods with their inherent advantages 
discussed in the section, Conjugate Direction Methods, can be used. Conjugate 
direction algorithms for computational solution of optimal control problems by 
Approximation in U have been formulated by Lasdon, Mitter, and Waren [29]; 
Sinnott and Luenberger [30], and Lynch [31]. In the section, A Conjugate 
Direction Method for Approximation in X 5 '” , we present a new technique for 
Approximation in U based on concepts from Hamilton-Jacobi theory and 
asymptotic stability. As an important special case, the state regulator problem 
is treated in detail. 

METHODS BASED ON THE SECOND VARIATION 

The following is a sketch of a second variation method that leads to an 
algorithm for Approximation in U . Because of the mathematical complexity 
of the results, we shall restrict this discussion to the special case of t^ 
specified and unconstrained terminal state. Details of the derivation and 
more general cases are given in the references [32,33]. 

Essentially, the method consists of minimizing a second-order approxi- 
mation to the Lagrangian equation (8) over the control space U and using a 
Riccati transformation to decouple the resulting boundary value problem. In 
fact, the entire development is similar to that given in the sections, Necessary 
Conditions, and State Regulator Problem in Chapter 2. Let us assume that 
the first-order necessary conditions, equations (35), (36), and (38), are 

(k) 

satisfied at the kth iteration by a control function u = u and corresponding 

state function x = x^ with increments <5u = 6u^ and 6x = 6x^ > 
respectively. Suppose the cost functional in equation (i) is written in the terms 
of the Hamiltonian functional, equation (34) , 



where the arguments of H and f have been omitted for brevity. Let J(u, x) 
be approximated by a constant term plus the sum of its first and second Frechet 
differentials (the first and second variations) at x and u . The procedure 
for obtaining the differentials is similar to that given in the section, Necessary 
Conditions, and the result is J(u, x) = J 0 + E(6u, 6x) where 

E ( 6u , 6x) = [^.SuJ+IJdx.H^dxJ+ljdx.H^du] 


+ 7 - fdu , H <5x1 +4- f 6u , H 6ul 
2 1 ’ux J 2 L uu J 


and J 0 is a constant. 

The increments 6x and <5u satisfy the linearized dynamical system 

<5x = f 6x + f 6u , <5x(t 0 ) = B , (43) 

xu u 

and superscripts denoting iteration are omitted. For convenience, we are 
using the Hilbert space inner product notation as defined in the Terminology 
section, Chapter 1. 

The design of a second variation algorithm is based on finding a control 
increment 6u which minimizes E(6u, 6x) subject to equation (43) . This 
"accessory minimization problem” is solved by the usual technique — equating 
the first Frechet differential E at <5u in an arbitrary direction £ to zero, 
integrating by parts, and applying the Euler-Lagrange Lemma (see the section, 
Necessary Conditions, Chapter 2). The result is as follows: 4 

8x = H Sx + H <5u , 6x(t 0 ) = B (44) 

Ax Au u 

6A = -H Sx - H 6u , <5A(ti) = B (45) 

xx xu 


4. Recall that H = f . 

A 
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( 46 ) 


0 = H + H <5u + H 5x + H SA 
u uu ux uA 

where SA is the adjoint function increment. 

Moreover, sufficient conditions for a weak minimum are the strong 
Legendre and the Jacobi (conjugate point) conditions: 


(a) 

H >0 
uu 


(47a) 

(b) 

H - H H" 1 H ^ 0 
XX XU uu ux 

. 

(47b) 


(An alternate approach is taken in the section, A Second Variation Method for 
Approximation in X , Chapter 4, where the minimization is performed in the 
state space X after applying an "inverse mapping.") 

If H is positive definite on T , its inverse exists, and equation (46) 
may be solved for the control increment as 

<5u = -H _1 (h + H 5x + H SA^ (48) 

uu V u ux uA / 

By eliminating <5u in equations (44) and (45) , we obtain the following linear 
boundary value problem in the product space X x X* x T : 

<5x = A 1 5x + A 2 <5A + v , Sx(t 0 ) = 9 (49) 

SA = -A s 6x - A t SA +w , SA(tj) = 9 (50) 


where 


A 1 = 


f - 

f H- 1 H 

A, = f 

H- 1 

f 

X 

u uu ux 

* u 

uu 

u 

H 

- H H- 1 H 

, V = f 

H 

H 

XX 

xu uu ux 

u 

uu 

i 


w = H H' 1 H 
xu uu u 
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In order to decouple equations (49) and (50), we introduce the "nonhomogeneous" 
Riccati transformation 

6A = P6x + q . (51) 

Then the solution to the boundary value problem may be expressed in terms of a 
matrix Riccati equation 

’ * 

-P = PA 1 + A 1 P-PA 2 P + A 3 , P(tj) = © , (52) 

and the linear equation 


-q = (PA 2 + A] )q + Pv + w 


q(tj) = 6 


(53) 


The sufficiency conditions [eq. (47)] insure that H -1 exists and that equation 
(52) possesses a solution on T . 

The basic second variation algorithm proceeds as follows: 

1. Choose a nominal control function u^ ; let k = 0 . 

2. Integrate equation (35) forward from x(to) and store the function 


x 


(k) 


(k) 


3. Integrate equation (36) backward from A(t-t) and store the function 


4. Compute the partial derivative matrices; check equation (47a) , 
compute IT^ , and check equation (47b) . 

5. Integrate equations (52) and (53) backward from P(tj) , q(t x ) , 

(k) (k) 

and store the functions P v , q v 

6. Compute the gain matrices 

Y = - H -1 (H +H P) 
uu ux uA 


z = -H- 1 (H + H a) 
uu u uA 
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7. 


<5x 


. 7. Integrate 6x = (f + f Y) <5x + z from 6x(t n ) = 6 and store 

(k) xu “ 

8. Compute the control update u ^ + *) = + Y<5x^ + z . 


(k) 

9. Let k — k + 1 , and repeate steps 2 through 8 until ||H' || 2 

is less than a predetermined positive number. 

Example 4 

Use of the second variation algorithm is now demonstrated by the 
computational example introduced in Example 2. The first iteration proceeds 
as follows: 

1. Let u ( °> - 0 . 

2. Integrate equation (35') forward from x(0) = i : 


(0) -t 
x = e 


3. Integrate equation (36') backward from A(l) = 0 


. (0) 

1 ( t-1 



A 

II 

to 1 
CD 

-e ) 

• 

4. 

H (0) 

= te -t 

H ( °> 

- H ( ? 

u 

Au 

uA 

H<°> 

XX 

= l 

H (0) 

H Ax 

= -1 

H (0 > 

II 

o 

II 

0 f = 

-1 

XU 

ux 

X 


H ( °> 

uu 

= i 

f = 
u 

i 
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5. Integrate equation (51) backward: 


P - 2P - P 2 + 1 = 0 , P(l) = 0 

for which a solution by separation of variables is 

p = (Nf2- 1) -Ml -~/2) e 2 ^" 


(Note: The function q is zero because v and w are zero. ) 

6. Compute the gain matrices Y = - P and z - -te 


-t 


Equation (43) may be integrated to obtain 6x and, hence, the control incre- 
ment 6u as specified in steps 7 and 8. 

7. Compute <5u^ , and update u^ : 


u 


( 1 ) 


' t-i -(t-l) . 

+e P(t) 


) 


Steps 2 through 7 are repeated until 


II H (k) II 2 = /(u (k) + A (k) ) 2 dt < 

U /-V 


The second variation algorithm requires integration of — n(n+9) equations 

& 

at each control function update and consequently suffers from the "curse of 
dimensionality. ” In other words, high-order problems may require an exces- 
sive amount of computation. In the last example, n = 1 , and five equations 
must be solved. Another disadvantage is the complexity of the computer pro- 
gram since the algorithm requires a significantly greater number of instructions 
than first-variation methods. In compensation for the latter difficulties, the 
second variation algorithm results in rapid convergence to a minimum. However, 
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the occurrence of conjugate points [violation of the Jacobi condition, equation 
(47b)] would invalidate the procedure, whereas first-variation methods could 
possibly be used to obtain a solution. Usually the full increment is not used, 
and adjustment of step length should be provided to insure stable descent. For 


example, let u 


(k+1) 


- u 


(k) 


a < k W k) 


where a ^ s 1 [32], 


A Conjugate Direction Method for Approximation in X' 

A detailed formulation of the boundary condition iteration algorithm using 
conjugate direction search in X* is developed for the optimal control problem 
of the section. Definition of the Optimal Control Problem, in Chapter 2. In 
contrast to former methods for boundary condition iteration involving the 
Newton-Raphson method [26,24], the inverse Jacobian matrix is obtained by 
iteration, and only first-order partial derivatives of the cost functional are 
required. This formulation has been reported by Gruver [9]. 


The most difficult aspect of setting up a boundary condition iteration 
algorithm is choosing a function of the terminal boundary conditions to be 
minimized. In the section, Descent Algorithms for Hilbert Space, we employed 
a positive definite functional involving the terminal adjoint function. A 
generalization of the latter to the case of t t free and terminal equality con- 
straints is to choose the Euclidean norm of an IR n+ ^ + ^ -valued function of 
terminal boundary conditions, equations (38) through (40). Suppose we define 

P 

the mapping E:XxX' xIR x T — IR by 


E(x , A , v , t) 


1 

2 


II cu l | 2 


l 


(54) 


where co 


is 


the IR 


n+p+i 


-valued function, 


co = ^H(x , A , u(x , A) , t) + £i^(x > t) , rj , A(t) - , t) v 

ipix , t) 

n_LryL i 

and || • || is the usual Euclidean norm on IR ^ . As before in the develop- 

l 

ment of the section, Descent Algorithms for Hilbert Space, d is the unknown 
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initial condition for the adjoint function in A in X* , v is the unknown con- 

* 

stant adjoint function with values in (EEr ) , and is the unknown terminal 

time. It is assumed that there exists a solution to equation (35) through (40). 
Then, if equations (35) through (37) are satisfied, a necessary condition for a 
weak minimum of J over the product space X x U x T is E(x,A,y,t 1 ) = 0 . 
We shall consider the simpler problem of choosing a minimizing sequence 




such that 



and perform the minimization by conjugate direction search using the quasi- 
Newton method of Fletcher and Powell [19], 

Next, the gradient is computed in a manner similar to that described 
in Example 2 in the section, Descent Algorithms for Hilbert Space. The 

gradient of E is the nt n+rH " 1 -valued function 

g(d, v, t) •= (E d , E v , E t ) | . (56) 

h 

Calculation of E and E may be performed explicitly. The ith element 

V i> 

of the IR -valued function E , is 

d 


n+p+1 n 


9 CO. 02 


9c o. 9 A, 


0x k 9d. 97^ 9d 


where x = (x ls . . . .x^) , A = (A lr . . . , A^) , and co = (co l5 . . . , <^ n+p+1 ) • 

Thus, evaluation of E, requires the n x n matrix functions x. and A, . 

’ d d d 

By differentiating equations (41) and (42) with respect to the initial condition 

d , we obtain the matrix "sensitivity equations," 

x , = H (x , A , u(x , A) , t) x, + H (x , A , u(x , A) , t) A, , 
d Ax d AA d 


x ( t 0 ) = © 


(58) 
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( 59 ) 


A , = -H (x , X , u(x , A) , t) x , - H , (x , A , u(x , A) , t) X. 
a xx a xA a 


A d (t 0 ) = I 


Equations (58) and (59) provide an accurate means for computing and, 

hence, the gradient of E in terms of (d , v , tj) . Knowledge of the gradient 
enables us to implement the conjugate direction methods described in the 
section, Descent Algorithms for Euclidean Space. 

The complete conjugate direction algorithm for Approximation in X* 
is as follows: Assume that equation (37) has been used to eliminate the depen- 
dence of equations (35) and (36) on the control function u as in equations (41) 
and (42) . Then 


Choose an initial it 


(0 ’ 0) =(d‘°>,,< 0) ,t 1 (0 >) ; leti = J.o. 

2. Integrate equations (41) and (42) forward from x(t 0 ) , A ( t^) = d^ 
until t = ti^ . Compute u(x , A) and store x^ (tj^^ , A^^t^^ 

3. Integrate equations (58) and (59) forward from x (t 0 ) , A, (t„) until 
t = tj , and store the functions x , A \ 


4. Compute and store g^ir^ ’ ^ = ( E , 

„(i , j) _ /^(i) ,.(i) ti (i) ^ _ V 


where it 

5. Compute it 


E , E 

v 


= Id , v 

(i , j) 


.) 


by the quasi-Newton method with Fletcher— Powell 


update. 


a. Choose = I . 

b. Compute s ( ^ = -H ( ^ g 


minimized. 


c. Determine a 

j (i , 3+1) 
a . it = 


(j) 


such that E 


(■“ ■“)■ 

(,« ■ » 


(j) (j) 


+ cx s ^ is 


j) +0! (j) B «> 


}. Compute g^7T (l 5 ^ J = ^E d , E^ , t^ 
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f. 


(j) 


2 


H 


( j +1) = H (j) + a (j) U-5- 


Ag (j) H (j) J[2 


[ 


Ag<»,s<»] [ 


Ag (j) , H (j) Ag (j) ] 


where 


Ag (J > = g( lr ,i ' i+1) ) -g^ 1 ’") . 


g. Let j — j+1 and repeat steps b through f until || g 
is less than a predetermined positive number. 


6. Let i — i+1 and repeat steps 2 through 5 until E 
than a predetermined positive number. 




is less 


Determination of the parameter a ^ in step 5c may be accomplished by a 
one-dimensional search procedure [14, 15,29]. In step 5f the quasi-Newton 

update as given by Fletcher and Powell [19] is used to update , an n x n 
matrix of functions on T which converges to the inverse Jacobian matrix of 

g^7r^ ’ ^ as j — <» . Notice that the indicated norms and inner products in 
step 5f are in the Hilbert space L^ + ^^(T) . 

The conjugate direction algorithm just presented may be applied to a 
wide class of optimal control problems and inherits the advantages of both the 
boundary condition method (moderate amount of storage) and the quasi-Newton 
method (rapid and stable descent). The algorithm requires fewer instructions 
than second-variation methods. Since the algorithm is based upon satisfying 
only first-order necessary conditions for a minimum of J , it does not depend 
upon the sufficiency conditions [eqs. (47a) and (47b)] second variation methods. 
Two very important advantages of the algorithm are the stability and rapid 
convergence resulting from the Fletcher-Powell update, step 5f, of the quasi- 
Newton method. The Newton-Raphson method locates the minimum of a 

quadratic functional J : IR n -*• IR in only one iteration and could be used to 

compute 7r^ ’ ^ in step 5. However, the iteration may be unstable for 
certain initial elements. Steepest descent methods, although very stable, 
may never converge if the eigenvalues of the Jacobian matrix are far apart. 
Quasi-Newton methods require at most n iterations to minimize a quadratic 
functional. The Fletcher-Powell quasi-Newton method results in rapid and 
stable descent even for highly nonlinear terminal boundary conditions and is 
acknowledged as the most powerful minimization technique currently available. 
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The most difficult aspect of using the algorithm is choosing the one- 
dimensional search procedure in step 5c. Quasi-Newton methods usually 
require accurate determination of the one-dimensional minimum. However, 
since each evaluation of E in step 5c requires forward integration of equations 
(41) and (42), direct search or higher techniques such as cubic interpolation 
may require too many evaluations of the functional E and, consequently, an 
excessive amount of computer time. It is suggested that a quadratic inter- 
polation scheme [ 1 5 J or golden section minimization [14] be employed, 
although the final choice of the best one-dimensional search must be resolved 
by numerical testing of the algorithm with nonlinear problems. 

EXAMPLE 5 


Use of the conjugate direction algorithm for Approximation in X “ 
will be demonstrated by the problem given in Example 2. The first iteration 
proceeds as follows: 

1 through 4 are identical to Example 2 with 7 ^ = d^ = 0 . 
5. Update 7r^ * ^ by the Fletcher-Powell method: 


a. H' 


b. s 


= - g (* (0 • °>) 


i / -41 1 4 2 t 

= - — 1 e - e 


-41 t 41 1 

0. 586 e + 2.414 e 


,(o) , (o) (o y 


c. Determine a such that Eyd + a s J is 
minimized (see Example 2, step 4, in the section, Descent Algorithms for 
Hilbert Space) . 

d. = *< 0 ’ 0 > +a < 0) 5 (0> = a ,0) s< 0) . 

e. Compute g^7r^ ’ from equations (58) and (59) as in 

step 4. 
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Steps 5b through 5f of the algorithm are repeated until 

n (g ").n* = f g(- (1 - J) ) 

{, . . , 0 

’ ,3 = 1 , 2 ,... 

/ (i oo) \ 

El 7r ’ 1 < a 5 . Notice that the first iteration of the quasi-Newton method 

in step 5, for = I , results in minimizing in the direction of steepest 

descent. However, successive step directions s^ are modified in step 5b 
by the matrix H ^ and eventually result in a Newton-Raphson iteration. 


dt < a 4 . After convergence of the 
, steps 2 through 5 are repeated until 


A Method for Monotone Approximation in U 

In this section, a sequence of functional approximations to the optimal 
cost functional is obtained by using concepts from the Hamilton- Jacobi theory 
to generate a minimizing sequence for J in the control space U . This 
technique was first reported by Gruver [13] and Puri and Gruver [34]. 5 
Similar results for a more general case were given by Leake and Liu [35] . 

An alternate proof for the case of the state regulator problem has been given 
by Kleinman and Athans [36] and Wonham [37], 

Consider the optimal control problem as stated in the section, Definition 
of the Optimal Control Problem, in Chapter 2. Let us assume that the terminal 
time is specified, and the terminal state is unconstrained. The following method 
is based on obtaining a (k+l)st control function approximation to the optimal 
control by minimizing a certain Hamiltonian functional containing the kth 
control function approximation. Let the kth approximation to the minimum 
(k) 

cost functional V ;XxUxT-*E be defined by 


V 


(k) 


t (lr ) 

J L (x , u , s) ds 
t 


(60) 


and the Hamiltonian functional H : X x U x X'" x T -* IR by 


H (x , u , , t) = L(x , u , t) + f V (k ' , f(x , u , t)1 

X LX J 


(k) 


( 61 ) 


5. This section and the next one contain an expansion and correction to 
References 13 and 34. 
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(k) 

By defining the Hamiltonian as in equation (6i) , we are treating V as an 

X 

adjoint function in X' 1 ' . This choice provides a convenient method for deter- 
mining the adjoint function from the cost functional. For example, suppose 

(k) 

that L is a quadratic functional on X x U x T . Then V has the repre- 
sentation = ifx, x j for some n x n matrix whose 

2L *1 (k) (k) 

elements are functions on T and, therefore, V' ' = P v 'x . Consequently, 

(k) X 

P' ' specifies a Riccati transformation as in equation (22) and, it can be shown, 

also satisfies a matrix Riccati equation. Further results on the correspon- 

( k ) 

dence of V v and the adjoint functions in X ‘ are given by Kalman [11] 

X 

although, for our purpose in the following section, the latter discussion pro- 
vides adequate motivation for use in the method that follows. 

Given the initial approximation u^ , suppose that we have obtained 

(k) ( k ) 

the kth approximation u v ' and the corresponding cost functional V v 

By the min-H condition a better approximation u^ k+ ^ to the optimal control 

u is obtained by minimizing H^x, u, V^ k ^ > overall ueU . In fact, the 

Maximum Principle [38] insures that u^ k+ ^ f s a better approximation even 
if the control function is restricted to a closed proper subset of U . In the 
latter case, however, there are no simple methods for performing the mini- 
mization since the gradient of H may not exist. 

We shall now show that the latter procedure for obtaining u^ + '*'^ from 
u v ' results in a sequence of cost functional approximations, which is mono- 
tone decreasing to the minimum cost functional. By definition of u^ +1 ^ we 
have 


„ / (k+1) (k) . , (k) 

H (x , u , V , t) = inf H (x , u , V , t) 

X T T X 

ueU 


which implies 


H (x , u (k+1) , V (k) , t) s H (x , u (k) , V (k) , t) 

X X 


( 62 ) 
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( k) 

Furthermore, it may be shown [11] that V v satisfies the "Hamilton- 
Jacobi" partial differential equation, 


H (x, u 


(k) „(k) 


t) + V 


(k) 


= 0 


( 63 ) 


Consider the control function u' and the induced state function from 

equation (2). Then by using equations (60) through (63) we obtain 


d_ (k+i) 
dt 


L (x , u (k+1) , t) 


j . (k+i) .. / . (k+i) (k) .. 

L (x , u , t) - H (x , u , V , t) + V 

V x 




tlx j] 


_d (k) 
dt 


(64) 


The second line in equation (64) follows since equations (62) and (63) imply 
that 


H (x , u 


(k+i) TT (k) 


t) + V, 


(k) 


0 


We are using the fact that the third line of equation (64) is the total derivative 
(k) 

of V ' ' with respect to t around the state function x as determined by 

u (k+1) . By integrating the first and last terms in equation (64) on (t, tj 

(k+1) (k) 

and using the fact that V 1 / (t[) = V v (t t ) =0, we obtain the inequality 


V (k+1 > s V ,k) 


k = 1, 2, . . . 


( 65 ) 


which proves that 



is monotone decreasing. 
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(o°) 

Since V also satisfies equations (62) and (63) and, therefore, for 


k = 1 , 2 , ... 


h(x , u , V^°° ) , t) + V t (00) ^ H (x , u , V^ k) , t) + V t ( 


(k) 


( 66 ) 


the limit V 


(o°) 


of the sequence 


1 v«) 


can be established by the same 


reasoning leading to equation (65) . It is more difficult to show, in general, 


that V 


(°°) 


A 


= V , the minimum cost functional. Leak and Liu [35] proved it 

(k) 

by assuming continuity of the operator that transforms the functional V' 

into . In the following section, we consider the special case of linear 

dynamical constraints and quadratic cost functional (the state regulator pro- 
blem described in the section, State Regulator Problem, Chapter 2). It is 
(k) 

shown that V ' may be expressed in terms of a positive definite matrix 


( k) 

P v , which satisfies a sequence of linear differential equations that converges 
as k — °° to the matrix Riccati equation (23). Since it has been shown that 
equation (23) possesses a unique positive definite solution [2], existence of 

the optimal control for the original problem implies that j j converges 

to the minimum cost functional. In the more general case, however, explicit 

(k) a 

conditions which insure that lim V = V are not yet available. 

k— oo 


At this point we shall assume that L(x, u, t) is positive semidefinite 6 . 
In addition, suppose tj — + ®° . Then we may guarantee the stability of 

equation (2) resulting from the approximation sequence J u^ J as follows. 

( 1 ) 

Assume that there exists an initial control function u eU such that 
equation (2) is uniformly asymptotically stable [39, p. 395]. Then it is 

shown in Reference 2 that is a Lyapunov functional for equation (2) 

with u = u (1) . Suppose that we have obtained the kth approximation u^ , 

(k) 

and the corresponding cost functional V ' is a Lyapunov functional for 

equation (2). Then by applying the previous method of choosing u^ k+ *^ by 

(k) (k+1) 

minimizing H(x, u, V' ' , t) , we can show that u v is also a "stable" 

X 


6. A functional L:XxUxT — IRis positive semidefinite if L(x,u, t) ^ 0 
for all xeX , ueU , teT with equality if, and only if, x = Q , u = 6 . 
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control. By assumption L(x, u 


(k+1) 


, t) is positive semidefinite, which 


implies that V^ k+1 ^ is positive definite. Moreover, -j|-V^ k+i ^ = -L(x, u^ k+: ^ , t) , 

which is negative semidefinite. Thus V^ k+1 ^ is a Lyapunov functional for 

equation (2) with u = u^ k+ *^ , which implies that equation (2) is uniformly 
asymptotically stable. 

The preceding development may be summarized by the following: 

Theorem : Given the cost functional [eq. (1)] and dynamical system constraints 

(k) 

[eq. (2)], let us define the mappings V :XxUxT-lby 
V (k ^ = j L^x , u (k ^ , sjds and H : X x U x X* x T — IR by 

H(x, u, , t) = L(x, u, t) + , f(x, u, t)l . Then given an initial 

(1) X * x !*i (k) \ 

u' , a sequence of control function approximations j u' J may be generated 


by minimizing H(x, u, V 


(k-1) 


x 


,t) , k = 2, 3, ... over all ueU, and this 


sequence possesses the following property of monotone approximation: 

a (°o) (k+1) (k) 

V <y ' <y l ’ < V v ' 

(oo) (k) A 

for k= 1, 2, ... where V = lim V exists, and V is the minimum 

k— *• °° 

cost functional corresponding to the optimal control function. Moreover, if 

L(x, u, t) is positive definite, tj — °° , and if there exists an initial u ^ 
such that equation (2) is uniformly asymptotically stable, the dynamical 

( k ) 

system, equation (2) is uniformly asymptotically stable for u = u , 
k = 2, 3, , . . 

The State Regulator Problem Revisited, I 

As an application of the method for monotone Approximation in U , 
let us consider a special case of linear dynamical system constraints 


x = Ax + Bu 


x(t 0 ) = c 


(67) 


and the quadratic cost functional 
J(u , x) 


= 1 ( [x . Ox] + [u . Ru] ) 


( 68 ) 
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I 


where A, B, Q, R are matrices whose elements are functions on T and 
satisfy conditions given in the section, State Regulator Problem, Chapter 2. 

Suppose that u^ is an initial control function. Then by the previous 

theorem, a minimizing sequence for J can be obtained by selecting u = u 
in U such that the Hamilton functional 


H ( x , u 



!([*• H + [" ’ R “] ,) + [<“' • Ax + Bu] _ 

(69) 


is minimized. Since U is an open set, a weak minimum of H can be 
found by equating its gradient with respect to u with zero. Hence, 


e = h ( x , u , v (k) , t) 

U X 


= B*V (k > + Ru (k+1) 

X 


(70) 


and since R > 0 we may solve equation (70) for u^ k+1 ^ = -R . To 

(k) X 

evaluate V v ' in terms of known quantities, we shall use an approach similar 

X 

to that described in Reference 2, which originally motivated the definition 

of . Assume that can be represented as = Trfx, P^xl 

(k) 21 1 

where P v ' is a positive definite, symmetric, n x n matrix whose elements 

are functions on T . Then the kth approximation to the optimal control 

function is 


(k+1) 

u 


-R- 1 


B 


O' 



k = 1 , 2, 3, . . . 


(71) 


Using equation (68) and changing the lower limit of integration to te[t 0 , tj , 
may be written as 





* (k-1) 

P x , 


-B* P (k_1, x 


1 ) 


ds 
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Differentiating the last equation with respect to t and substituting 

. ji ; / 1C -- 1 3 

A - BR -1 B ' P / x for x , we obtain the following linear matrix differential 
equation: 

-P< k > = P^ k L^ + (a^)* P (k > + Q (k > , P <k) (tj) = • 

(72) 

where 


A (k) = A - BR- 1 B*P (k 1} 


=Q+p^ k ^ BR- 1 B*P^ k ^ 


(73) 

(74) 


Equations (71) through (74) define a recursive solution for the optimal control 
functions which may be summarized in the following algorithm: 

1. Choose an initial control function u^ and determine P^ ; 
let k = 1 

(k) (k) 

2. Compute A v ' , Q v from equations (73) and (74). 

„ _ , . * (k) fk) (k) . . (k). *_(k) (k) . , 

3. Integrate -P 7 = P A ' + (A ) P + Q backward 

(k) (k) 

from P v ' (tj) , and store the function P v 

(k) 

4. Integrate x = A' 'x forward from x(t 0 ) , and store the 
function x^ . 

5. Compute the next control function according to 
u (k+1 > - -R-‘B*P< k) X (k) . 

6. Let k -** k+1 and repeat steps 2 through 4 until 

|| ;g : ‘ c y( k ) + Bu^ k+1 ^ || 2 is less than a predetermined positive number. 

X 

In general, step 3 involves integrating a time-varying linear differential equation. 
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If a closed-loop control function is desired, step 4 may be replaced by com- 
putation of the feedback gain matrix R -1 B*P^ . For the systems of order 

n 2 3 , finding the required initial matrix P^ is nontrivial. The pole- 
shifting technique used in References 40 and 41 is an effective method for 

constant coefficient systems. However, a general method for selecting P^ 
is not available at the present time. 

The greatest advantage of the algorithm just presented is that equation 
(72) specifies a linear matrix differential equation that must be solved to 

compute the new control function u^ + ^ . Methods for solving equation (72) 
are given in Appendix A. (In contrast, the standard approach given in the 
section, State Regulator Problem, Chapter 2, involves solving a nonlinear 
matrix Riccati equation for which a solution may not even exist. ) The 
algorithm requires a modest amount of computer instructions and computation 

time. It was recently shown by Kleinman and Athans [36] that | j also 

possesses quadratic convergence, whereas most other algorithms for computing 
the optimal control, such as ASP or Runge-Kutta integration of equation (72) , 
display only linear convergence to the minimum cost functional. Although 
further testing of the algorithm with high-order linear time-varying systems 
having finite time interval (t t < °°) is needed, numerical studies by this 
author and application of the algorithm by others [36,41,40] have indicated 
that the method for monotone Approximation in U described in this section 
is the most efficient technique for the case of constant coefficient systems 
and infinite time interval. 


EXAMPLE 6 

Given the quadratic cost functional 
1 *1 

J(u,x) = — j (qx 2 + ru 2 ) dt 
2 0 


(75) 


and the linear first-order system 

x = -ax + u , x(0) = c , (76) 

where q s 0 , r > 0 , and < °° . The problem is to determine a minimizing 
sequence for J . Computation of an approximation to the closed-loop optimal 
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control function by the previous algorithms for monotone Approximation in U 
proceeds for the first iteration as follows: 

1. Choose u (1) = 0 , which implies P^ = 0 . 

2. Compute = -a , = q . 

3. Integrate equation (72): 


-P* 1 ^ = -2aP*^ + q 


P (1) (ti) = 0 


which has the solution, 


P^ = ^ (i - exp (-2a (t t - t) ) ) 


4. Update the control function: 


u< 2 ’ , ipl'IW 

r 


= ( i - exp(-2a (tj - t))) x 


(i) 


The second iteration requires integrating the time-varying linear equation 


P l2 > ■= -2A ^ P^ - Q ^ 


p (2, (t,) = 0 


where 


A^ = a +~ — ( 1 - exp (-2a (tj - t)) ) 


2ar 


= q + ( 1 - 2 exp ( - 2a (tj - t) ) + exp(-4a (tj - t)) ) 


4a r 
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The solution to the last equation is 


p(2) _ ^ ex p(_2an (t t - t) ) 

n=0 n 


where 


p(2) = q( 4a 2 r + q) 

0 4a(q + 2a 2 r) 

p(2) . 

1 2a ( 4a i! r + q) 

(2) = q ( 4aP l 2>+q ) 

2 4a( 6a 2 r + q) 


and 


oo 


l 

n=0 



0 


Therefore, the resulting approximation to the optimal control function 
is 

o° 

u (3) = ~~ Tj p{ n 2) exp(-2an(t 1 - t) > x (2) 
n=0 

Suppose that tj — + <*> . Then if a > 0 , a minimizing sequence of asymptoti- 
cally stable functions for J is as follows: 


u 


( 1 ) 


0 


(77) 
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(78) 


u 


( 2 ) 


2a 


x 


( 1 ) 


(3) _ 4a 2 + 1 (2) 

U 4a (2a 2 + i) X 


EXAMPLE 7 

Given the cost functional 


J (u , x) = “ f (xf + u 2 ) dt 
2 0 

and the linear system 

xi = x 2 , x^O) = CJ 

x 2 = -X!-ax 2 +u , x 2 (0) = c 2 

where a > 0 , we apply the algorithm of this section to find a minimizing 
sequence for J . Let us define the state function x = (x 1} x 2 ) and the 
matrices 



Q = diag(l ,0) R = 1 

A minimizing sequence u^ is obtained from 

/ k) 

The constant matrix P J is the solution of the linear matrix equation 
P (k >A (k) + (A (k >)"p (k) + Q (k ^ = © , 


(79) 


(80) 

(81) 

(82) 


(83) 

(84) 
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where 


A^ = A - BR -1 B*P 


= Q + P (k " 1) BR -1 B*P (k-1) 


Since a > 0 we may choose u^ =0 ; hence = © , and from 

equations (83) and (84) the result of the first three iterations is 


u 


u 


u 


( 1 ) 

( 2 ) 

(3) 


x * 


10a 2 + 3 


12 Xl - 1 2a ( 2a 2 + 1 ) Xz 


EXAMPLE 8 


Given the cost functional 


(85) 

( 86 ) 

(87) 


J(u , x) = ^ f (qn x? + q 22 x^ + ru 2 ) dt 
2 0 


( 88 ) 


where q^ ^ 0 , q 22 ^ 0 , r > 0 , and the linear third-order system 


Xl 

= x 2 

(89) 

*2 

= *3 

(90) 

*3 

= -0. lxj - 1 . 2x 2 - 2. lx 3 + u 

(91) 
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determine a minimizing sequence for J . Applying the previous algorithm, 
we define the state function x(t)eIR 3 , control function u(t)eIR , and the 
matrices 



Q = diag(q tl , q 22 > 0) R = r 

The problem was programmed for the digital computer, and representative 
results are shown below. In each case the closed-loop approximation to the 


optimal control is 

u = 

3 

- Yj k x * 

n=l n 




on 

0.22 

r 

ki 

k 2 

k 3 

min 

cost 

l 

0 

0. 5 

1. 3177 

1. 8090 

0. 7334 

2. 0730 

i 

0 

1 

0. 9050 

1. 3215 

0. 5558 

2. 4141 

i 

0 

2 

0. 6141 

0. 9461 

0. 4104 

2. 8252 

i 

1 

0. 5 

1. 3177 

2. 1787 

0. 8610 

2. 3351 

i 

1 

2 

0. 6141 

1. 0756 

0. 4615 

3. 0103 


In each case, four iterations were required to obtain convergence of the 
algorithm. As a check on the computer results, the exact solution was 
obtained for q^ = 1, q 22 = 0, r = 1 using a spectral factorization method 
devised by Puri and described by Puri and Gruver [34]. For this case the 
exact solution to three-decimal accuracy is: 

kj = 0. 905 k 2 = 1. 32 k 3 = 0. 556 
which compares favorably with the approximation shown above. 


52 



Further examples of the latter approximation method to the state reg- 
ulator problem have been given recently by Me Lane [41] and Yu , Vongsuriya, 
and Wedman [40] . In the latter reference, the authors treat the optimization 
of a power system involving an eighth-order constant coefficient system with 
two control inputs. Analog simulation was used to check the computer solution, 
and the results appear favorable. 


CHAPTER 4 


APPROXIMATION IN X 


In this chapter, we formulate a new class of computational methods 
for the solution of the optimal control problem stated in the section, Definition 
of the Optimal Control Problem, Chapter 2. Given the cost functional [eq. (1)] 
subject to dynamical system constraints [eq. (2) ] , we seek a minimizing 
sequence for J by Approximation in X . Central to the theory is a represen- 
tation for the mapping from X into U . A new second variation algorithm is 
developed for the nonlinear control problem. In contrast to previous second 
variation methods, we avoid the use of adjoint functions and the Hamiltonian 
formulation. At each control function update, the control increment is obtained 
in terms of the state increment in a least-squares sense. An accessory 
boundary value problem is obtained, the solution of which specifies the control 
function update. As an important special case, the state regulator problem is 
treated, and several computational examples illustrate the use of the technique. 


The Inverse Mapping 

Let us assume that the dynamical system [eq. (2)] may be linearized 

(k) 

about an initial control function u = u' ' in U and corresponding state 
( k) 

function x = x' in X . Hence, 

<5x(t) = f (t)6x(t) +f (t) <5u(t) , <5x(t 0 ) =0 , (92) 

X u 

(k) (k) 

where 6x = 6x v ' in X and 6u = 6u v ' in U are increments at x and u, 
respectively. For brevity, the arguments x, u of f and superscripts 
denoting iteration are dropped. 

If the matrix f (t) were invertible and nonzero for all teT, we 

u -i 

could solve for 6u(t) simply by premultiplying equation ( 92) by f (t) . In 

general, however, f (t) is not invertible and the pseudo inverse [42], 

f > K 0 ‘> v ”)" 1 f ; ;<•> • 
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yields the least-squares solution for Su(t) [in terms of Sx(t) ] to the 
"overdetermined" system [eq. (92) ] , which has more equations than unknowns. 
If f (t) is invertible, the psuedo inverse reduces to the usual inverse. 

In general, f^(t) may be zero for some te T. If f *(t) f (t) ^ 0, it is 

convenient to define the e-pseudo inverse of f (t) as 

u 

g e (t) = ( f u (t) f u (t) + eI )~ lf u (t) ’ < 93 > 

where e is a small positive (real) number. The parameter e insures 
that the indicated inverse in equation (93) exists. If f* (t) f^ ( t) > 0 for all 

teT, we shall choose e = 0. 

Based on the above remarks, equation (92) may be solved for the 
control increment in terms of the state increment as 


6u(t) = g £ (t) (£ 6x)(t) 


(94) 


where 


(£ <5x)(t) = <5x(t) - f (t) <5x(t) 

X 


(95) 


At this point of the development, we must restrict the function <5x to be at 
least continuous and possess a derivative in X . This condition is satisfied 
by assumption 1 of the section, Basic Assumptions, Chapter 2. Let us define 
the residual 

(p(6u))(t) = (£ <5x) (t) -f u (t) <5u(t) . 

The functional ||p(6u) |l 2 is a measure of the error between the acutal and 
desired solution to equation (92). For example, if f is invertible, then 

||p(6u)|| 2 = 0. It is possible for ||p(6u)|| 2 = 0 even if f is singular. 

An important example of the latter is a linear time-varying system having 
a single control function input. Use of the usual pseudo inverse (e = 0) in 
determining the control increment results in minimizing the functional 
||p(6u) || 2 . The following theorem specifies the function that is minimized 
by the use of the e-pseudo inverse: 
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Theorem : Given equation (92), suppose that the function £dx and the matrix 

operator f are specified in advance, where £ dx e X , due U , and 

induces a bounded linear mapping from U to X . Let e be chosen such 
that det(f*f^ + el) ^ 0 for all teT and let !lp e (6u)|| 2 be defined by 

1 1 p^ ( du) 1 1 2 = 1 1 f^ du - £ dx 1 1 2 + e 1 1 du 1 1 2 . (96) 


Then, 


where 


p ( du) || 2 < || p (du) II 2 for all <5u e U , 
6 6 


du = (f*f + el)' 1 f* £ dx . 
u u u 


Proof: 1 1 p ( du) || 2 = II f du - £dx || 2 + e II du || 2 
e u 

= £f^ du - £6x, f^ du - £dxj+ e £du, duj 
= [f* f u du, du]+ [£dx, £dx J 
- £ f u du , £dx J - £ £dx , f^ du J -s- * £du , du J 

= [f>f u 6 u ,6„]- 2[du, ( rf u+ €l)6v] 

+ || £dx || 2 + e £du, du J 


where 


dv = (f*f + e I) _1 f* £dx 
u u u 
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Finally, by collecting terms, ||p e (<5u)|| 2 may be written as 

[| P e (6u) II 2 = £(f* f u + e I) (Su - 6v) , (<5u-Sv)J 

[ (f*f + el) <5v , <5v1 + 1 1 £ <5x 1 1 2 . 

u u J 


( 98 ) 


By assumption ||£6x|| 2 , a fixed number, is finite. Since 
f * f + e I > 0 , ||p e (du) || 2 is minimized by 6u = 6v . QED 

In summary, we have obtained in inverse mapping in the sense that 
the control increment may be expressed in terms of the state increment by a 
linear mapping of X into U as defined by equation ( 94) . Under certain 
conditions as specified in the last theorem, use of the e-pseudo inverse 
results in minimizing the square of the norm of the residual plus an additional 
term involving the control increment "energy." If f (t) =£ 0 for all t e T , 
we shall choose e = 0 . 


A Second Variation Method fo Approximation in X 

In the section, Descent Algorithms for Hilbert Space, Chapter 3, we 
described the well-known "direct" second variation algorithm based on 
minimizing, in the control space U , the sum of the first and second variations 
of the cost functional. In this section, we exploit the inverse mapping, defined 
in the preceding section, for the purpose of developing a new second variation 
algorithm that is based on performing the minimization in the state space X . 

Suppose the terminal time is specified, and the terminal state is 
unconstrained. Then let the cost functional 


J ( u , x) = 


f L (x , u , t) dt 
to 


(99) 


be approximated by a constant 
Frechet differentials of J at 
equation (1) with increments 


term J 0 plus the sum of the first and second 
u = u (k > and x = x^ , which satisfy 
= 6u^ and 6x = 6x^ , respectively, 
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which satisfy equation (92). This sum is computed by taking the following 
Gateaux differentials: 


E ( Su , 6x) = 


80?! 

J_ 

2 0Xj 3 x 2 


J L (x + oq 6x, u + oq 6u , t) dt 
t° 

0 2 / h 


a, = 0 


J L(x+ oqdx + 0*2 5 2 x, u + aqdu 
to 


+ aq <5 2 u , t) dt 


oq = 0 
a 2 = 0 


where 6 2 u and d 2 x are second-order increments at u and x , respectively. 
For brevity, arguments of L and superscripts denoting iteration are 
dropped. By expanding L in a second-order Taylor series and performing 
the indicated differentiation, E may be rewritten using the Hilbert space 
inner product notation as follows: 

E< 6u , Sx) =. , foj + [ L u • 6u ] + ~2 [ 6x ’ L xx 6X 1 

+ T [® I ’ L xu' 5u ] + ~H 5 u ' L ux' 5x ] “00) 

+ -4- f <5u , L dn 1 
2 [ ’ uu J 


The inverse mapping defined by equation (94) is now used to eliminate explicit 
dependence in equation (100) with respect to 6u . Thus we obtain 
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E(* , 6x) = J L x , 6x J + g e (dx - f x 6x J 

+ T I* 1 ' L xx 6,1 ] + T [ 5x ’ L xu g c (5i ' f x Sx) ] 

+ T [ g e < a4 - f x 6x) ’ L u x ,5x ] 

+ T [ g £ (6i - f x 5x) • L UU g s < ' 54 ■ f x 60 ] 

= [l x , &=]+ [ L u , g e ax] - [l u , g £ f x 6x] 

+ 4" Tfix, L <5x1+ 4- ffix, L g 6x 1 

2 L xx J 2 l xu 6 c J 

- — fdx, L g f 6x1+ 4" fg 6x, L 6x1 

2 L xu e x J 2 |, e ux J 

- 4" f g f 6x , L 6x 1 + 4" f g 6x , L g 6x 1 

2 l e x ’ ux J 2 [ °e uu s c J 

- 4" r S <5x, L g f 6x 1 - 4- Tg f <5x, L g 

2 I, & € uu & e x J 2 [ e x uu 6 

+ 4” Tg f 6x, L g f 6x1 
2 l e x uu e x J 


( 101 ) 



The development of the second variation method for Approximation in 
X is based on finding a state increment 6x in X which minimizes E . 
Since E inequation (101) is an explicit function 6x, this minimization may 
be performed without the use of adjoint functions and the Hamiltonian formu- 
lation used in the section, Descent Algorithms for Hilbert Space, Chapter 3. 
A necessary condition for a weak minimum of E over the state space X is 
that its first Frechet differential at 6x with increment £ vanish. Since it 
is assumed that the Frechet differential of E exists, it may be computed by 
using the Gateaux differential. The latter is linear in £ an d can be written 

in the form dE(6x;£) = £e £ J , where E^ is called the gradient of E 

with respect to 6x . Rather than compute the Gateaux differential as in the 
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section, Necessary Conditions, Chapter 2, we shall take the following 
alternative approach and (1) simply compute E^ ; 7 (2) use the definition 

of adjoint from the section, Terminology, Chapter i, which permits moving 
an operator from one side of an inner product to the other side; (3) integrate 
terms involving fix by parts, which results in "boundary terms" involving 
dxltj) and 6x(t 1 ) only, since <5x(t 0 ) = 0 . For example 


[V g e ai ]* [ g c L u’ 6i ] 




d_ , 
dt 

Then if we define 


_j , 

where 2D = I— is defined as the matrix differentiation operator on C X (T) 
t dt 


h = g L g , 
e uu e 


the differential of E is computed as follows: 

dE(6x;£)= f(L -^ (g*L ) -f*g*L 
lx t e u X e u 

+ L 6x + L g fix - 2Du L g )*fix) 
XX XU e t xu e 


- ( L g f + ( L g f ) * ) fix 
XU e x xu e x 


y) ( h fix) - — f * h fix 
t 2 x 


+ — 2D ( h f fix) - -4" f * h fix 
2 t x 2 x 

+ — (hf fix) + f* hf fix » £ 1 
2 t x x x J 

+ [< g c L u ) lt 1 +(g ‘ £ L xu- L xu g s ) lt, 5x(t ‘ ) 

+ (hf )| fix(ti) , |(t,)l 

X Itj Ji 

7. Note: The gradients of Jx, AyJ^ and Jy, AxJ with respect to x are 
Ay and A*y, respectively. 
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By collecting terms involving 6x and its derivatives, dE(<5x;|) may 
be written as 


dE( 6x ; {) = f L - & (g* L ) - f* g L 
lx t g u x g u 


+ (L (L g )* - (L g f + (L g f )*) 

xx t xu e xu g x xu g x 


-^(hf ) + f* hf )6x+ ( ( L g-(L g )*) 
t x x x xu g xu g 


-S>h. - f*h+ hf )5x 
t x x 


- h 5x , £ 1 + boundary terms 


( 102 ) 


Suppose we arbitrarily set the boundary terms to zero. Then by the Euler- 
Lagrange Lemma in the section, Necessary Conditions, dE(6x; £) = 6 implies 

that the state increment 6x = Sx ' satisfies the following linear second-order 
differential equation: 


h <5x = F <5x + G <5x + k 


where 8 


h = 


g* L g 
g uu e 


F = L g - g* L* + hf -f*h-^>h 
xu e e xu x x t 


G = L -L g f -f*g*L* + f* hf -3>.( g*L* + hf ) 
xx xu e x x e xu x x t e xu x 


k = L x- f i g * L u'^t (g J L u) 


(103) 

(104) 

(105) 

(106) 

(107) 


8. Recall that g = (f'" : 'f + el) -1 f * and 2D 

°e u u u t 


I 


_d_ 
dt * 
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and subject to the following set of "self-adjoint" boundary conditions: 


<5x(t 0 ) = 0 (108) 

h(ti) <5x(ti) = N(t t ) 6x(t!) + r(tj) , (109) 


where 


N 


hf - L 
x xu 


g + g * L 
e e xu 


(HO) 


r = g*L (111) 

e u 

( k ) 

To compute the control increment 6u = 6u ' from equation (94), we must 
solve equations (103) through (ill). Because of the similarity of the latter 
two-point boundary value problem to equations (49) and (50) of the section, 
Descent Algorithms for Hilbert Space, Chapter 3 (in which we obtained second- 
variation necessary conditions for Approximation in U ) , we shall refer to 
equations (103) through (111) as the "Accessory Problem." Computational 
solution of the Accessory Problem is discussed in Appendix B. The solution 
for the special case of the state regulator problem is given in the section, 

The State Regulator Problem Revisited, H. Sufficient conditions for the 
existence of a solution to the Accessory Problem, and hence existence of a 
weak minimum of E , follow from the Calculus of Variations and involve the 
classical strong Legendre and the Jacobi conditions f6,7]. 

EXAMPLE 9 

Consider the cost functional 


J(u, x) 



( 112 ) 


and the first-order linear dynamical system constraint 


KjZ 


x = u , x( 0) = 1 


(113) 



For this simple case 


L = u + 1 
x 

L = x + u 
u 


o 

II 

o 

11 

XX 

X 

L = L = i 

f - i 

XU XU 

u 


L = 1 
uu 


Since (f*f^) (t) > 0 , t<=[0 , 1] , we may choose e = 0 . Then g = i and 

h = 1 G = 0 N= 0 

F= 1 k=u+i-x-u r = x + u 

Therefore, equation (103) becomes 

Sx = u + 1 - x- u (103’) 

subject to equations (108) and (109), the "self adjoint" boundary conditions 

<5x( 0) = 0 (108') 

6x( 1) = u( 1) + x( 1) . (109*) 

A general algorithm for the computational solution of problems of 
this kind is given in the next section. 


The Basic Computational Algorithm 

The steps in computing optimal control functions by the methods 
described in the sections, The Inverse Mapping, and a Second Variation 
Method for Approximation in X , may be summarized in the following 
algorithm: 

1. Choose an initial control function u^ ; let k = 0. 
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2. Integrate x = f(x, u, t) forward from x(to) and store the 

- + . (k) .(k) 

functions x , x 

3. Compute the following partial-derivative matrices evaluated at 

x (k) and u (k) :f,f,L,L,L,L,L,L . 

X u u x’ ux XU UU XX 

4. Solve the Accessory Problem, equations (103} through (11 J), for 

6x^, (see Appendix B) . 


5. Compute the control function increment 

* (k) . _ .(k) (k) 

<5u = g (6x - f 6x ') . 

e x 

6. Update the control function according to 

(k+1 ) (k) (k). (k) , (k) . 

u = u + a ou , where a 5 1 isa constant. 


7. Let k-*k+l and repeat steps 2 through 6 until 

, (k) (k) . T/ (k+1) (k+1) , . . , . J 

|J(u , x )-J(u , X )| is less than a predetermined positive 

number. 


As was mentioned for the second variation algorithm described in the section, 

Descent Algorithms for Hilbert Space, Chapter 3, the full increment ou^ 

(k) 

should not be used, and the constant a K is included in step 6 to provide stable 
descent. 


EXAMPLE 10 

As an example of the second variation algorithm of this section 
consider the problem introduced in Example 9 of minimizing the cost 
functional 


J(u, 




dt 


subject to a first-order linear system 
x = u , x( 0) = 1 
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It was shown in the previous example that the Accessory Problem is 


<5x = u + 1 - u - x 


( 103 ’) 


subject to "self adjoint" boundary conditions 


<5x( 0) = 0 

6x( 1) = u( 1) + x( 1) 


(108') 

(109 1 ) 


The steps in computing the optimal control are as follows: 


1. Let u (0) = 0 


( 0 ) 


2. Integrating equation (113) yields x = 0 . 

1 ? 

3-4. Solve equation (103 1 ) by letting 6x = — V * + cqt + c 2 and 

5 

evaluating constants c* = - — , c 2 = 0 by using equations (108*) and (109'): 
(4. 17'-4. 18') : 




5. Compute the control increment from equation (94): 


<5u = t — 

4 


6. Update the control function: 


u 


(1) 


u 


( 0 ) 


+ 6u 


(0) 



In this example, J is a true quadratic, and the system is linear. Therefore, 

the differential is also linear in x and u , and the optimal control u = u^ 
has been obtained in a single iteration, as may be verified by comparison with 
the solution given in Reference 43. 
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EXAMPLE 11 

Next, we consider the minimization of the quadratic cost functional 
1 2 

J(u,x) = — / ( x+ “J u 2 ) dt 

subject to a first-order nonlinear system 
X = -X 2 + u , x( 1) = 1 

The Accessory Problem for this example is 

6x = ( 4x 2 - 2x) Sx + 1 + 2xu - u 
= ( 6x 2 - 2u) <5x + 1 + 2xu - u 

subject to boundary conditions 
<5x( 1) = 0 

6x(2) = -2x(2) 6x(2) + u(2) 

The steps in computing the first approximation to an optimal control are as 
follows: 

1. Let u<°> = 0 . 

2. Integrating x = -x 2 yields x^ = . 

3-4. Solve the Accessory Problem, 

6x = 6x + 1 

6x( 1) = 0 
6x(2) = - 6x(2) 


66 



The general solution to the latter equations may be obtained by integration 
and is follows: 


. 1 o 3 -2 

6x= lo* + ^ t 


i_ 

4 


t 2 


5. Compute the control increment from equation (94): 



6. Update the control function: 


u 


( 1 ) _.,(<» 


= u 


+ <5u 


( 0 ) 



In this example involving a nonlinear system, we obtained a second-order 
linear, time-varying equation for the Accessory Problem which could be 
solved analytically in closed form. In general, the accessory problem must 
be solved numerically on a digital computer using the methods given in 
Appendix B. 


The State Regulator Problem Revisited, 1 1 

Consider the state regulator problem of minimizing the quadratic 
cost functional 


J (u , x) = 


•r([ x ' Qx ] ' [ u ’ Ru ]) 


(114) 


subject to linear dynamical system constraints 


x = Ax + Bu , x(t 0 ) = c 


(115) 


where Q, R, A, B are constant matrices on T , Q s 0 , R > 0 , and the 
pair (A, B) is completely controllable. A computational solution to this 
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problem is now formulated by Approximation in X , using the methods described 
in the sections, The Inverse Mapping, and A Second Variation Method for 
Approximation in X . 

Since the cost functional is a true quadratic and the system is linear 
in x and u , the differential is also linear in x and u . Therefore, as 
shown below, equations (103) through (11) may be solved exactly in a single 
iteration. Consequently, we shall replace 6u and <5x by u and x , respec- 
tively, and derive the Accessory Problem for this special case. By means of 
the inverse mapping defined by equation (94) , the cost functional may be 
expressed explicitly in terms of the state function. Since B* B > 0 for all 
teT , the parameter e may be chosen zero. Then the cost functional may 
be written as 

J ( * , x) = ~ f[x , Qx] + [b + £x , RB + £x] \ (116) 


where 


B + = ( B*B) “* B* 

(117) 

£x = x - Ax ' . 

(118) 


A necessary condition for a weak minimum of J over the state space X 
is that its first Frechet differential at x with increment £ vanish. By 
performing the same steps as used in equations (103) through (111), we 
obtain the following linear second-order differential equation 

hx = (hA - A*h) x + (A*hA + Q) x , (119) 


where 


h = (B + )*RB + 

and subject to the "self adjoint" boundary conditions 


( 120 ) 


x(t 0 ) = c (121) 

h x(t t ) = hAx(tj) . (122) 
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After solving equations (119) through (122) for x and x , the optimal control 
function u may be recovered, using the inverse mapping, as 

u=B + (x-Ax) . (123) 


For the solution of the Accessory Problem, equations (119) through (122), it 
is convenient to consider the two cases: (1) h > 0 ; and (2) h S: 0 . When 

applicable, the solution for Case 1 is easier to implement. In the following, 
we shall obtain closed-loop control functions by using the Riccati transforma- 
tion. In Appendix B, both open and closed-loop control functions are obtained 
for the more general problem of the section, Definition of the Optimal Control 
Problem, Chapter 2. 

CASE 1: h > 0 

Since R > 0 , the condition h > 0 implies that the Euclidean 
dimension of the control space U and the state space X is equal; that is, 
dim U = dim X . Consider the Riccati transformation 


x = Pj x , (124) 

where Pj is an n x n matrix whose elements are functions on T . By 
substituting in equation (119) for x and x using equation (124), we find that 
P t satisfies the following matrix Riccati equation: 


-hPj = hP|+ (A*h - hA) Pi - (A*hA + Q) , P^tj) = A . (125) 

From equations (123), (124) and (117), the optimal control for this case is 

u= (B*B) _1 B*(Pj - A) x , (126) 

where P A is obtained by solving equation (125). 

CASE 2. h > 0 . 

The condition that h a 0 implies that dim U £ dim X . Let 
equations (119) through (122) be rewritten as the following pair of first-order 
vector differential equations: 

k a = A 1I x a +A 12 x b (127) 
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.b 

x 


(128) 


= A 21 x & + 


, 22 


where 


x — x — (xj , , xj , x x (^ n+ .j * • • * * X 2 n ^ 


n 


x. = x. , i = 1, 2, . . . , 2n - 1 
i+l i 


with boundary conditions 

x a (t 0 ) = c (129) 

C x a (tj) = x b (ti) . (130) 

Matrices A 1 " 1 , C are n x n with elements that are constant on T . Consider 
the Riccati transformation 

x b =P 2 x a , (131) 

where P 2 is n x n and has elements that are functions on T . By using 
equation (131) in equations (127) through (130) we obtain the following matrix 
Riccati equation: 

-P 2 = P 2 A u - A 22 P 2 + P 2 A 12 P 2 - A 21 , P 2 (t t ) = C . (132) 

Having obtained P 2 by solving equation (132) , the closed-loop optimal control 
for this case is obtained from equations (123), (131), and (117) as 


u= (B*B) " 1 B*(A 11 + A 12 P 2 -A) x 


(133) 
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EXAMPLE 12 


Consider the minimization of 


J(u , X) = f ^X 2 4 ““ u 2 jdt 

(134) 

subject to the first-order linear system 

x = - — x + u , x(0) = 1 

(135) 

Using the methods described for the state regulator problem of this section, 
the Accessory Problem is 

9 

X= — X 

(136) 

x(0) = 1 

(137) 

x( i) = — x( 1) 

(138) 


Since h = 1 is positive, Case 1 of this section applies. Solution of equations 
(136) through (138) by the Riccati transformation x = P(t)x results in the 
Riccati equation 


• 9 9 1 

- p =P ! - T , P(l) = --|- 


(139) 


The latter equation may be integrated backward by replacing the independent 

dP dP 

variable t by -s and letting — = - — r- . Then by separating variables 
and integrating, we obtain 


P= -f-coth (- -f-t+4,) 
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mi mini 


where 4 1 is adjusted such that P(l) = - ^ . Hence, from equation (126) , 
the optimal control is 

“ = ( P ""2') X '('^' + "| 00th ("4 t+4 ‘)) X ' <140) 

To check the validity of the solution just obtained, we employ the 
standard method for solving such problems as described in the section, State 
Regular Problem, Chapter 2. The Riccati equation (47) is 

p = p + P 2_ 2 , p(l) = 0 

which has a solution 

P = 2~ + ~2* co *k ( f” ‘ + 4 2 ) • 

where £ 2 is adjusted such that P(l) =0. Therefore, the optimal control 
is 



which is identical to equation (140). 
EXAMPLE 13 

Given the quadratic cost functional 


J(u, x) = — f (xl+ u 2 )dt 
^ 0 


and the second-order linear system 
Xj - x 2 , xj(0) = 1 


(141) 


(142) 
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x 2 = - x t - x 2 + u , x 2 (0) = 0 , 


(143) 


we seek the optimal control by Approximation in X . 

A necessary condition for a minimum of J in X is given by the 
system of equations (119), (121), and (122), which for this problem is 


0 = x 2 + 2xj + x 2 (144) 

x 2 = -Xj + Xj + x 2 , (145) 

where 

x 2 ( 0) = 0 (146) 

Xj ( 0) = 1 (147) 

x 2 ( 1) 4- x t ( 1) + x 2 ( 1) = 0 (148) 


For the preceding boundary value problem the matrix 


h = 


0 0 
0 1 , 


from equation (120) is positive semidei'inite. Thus, Case II of this section 
applies and equations (144) through (148) are written as 
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subject to the boundary conditions 


x a (0)=(i,0) (150) 

x b (i)=Cx a (i) , (151) 

where 



Equations (149) through (151) may be solved on the finite interval T by 
using the Riccati transformation (the sweep method) described in Appendix B. 
The optimal control is given by equation (133) , where elements of the 2x2 
matrix P 2 satisfy equation (B. 23) , which results in the following equations: 


-Pll = 

- P21 + Pll Pl2 > 

Pud) = -i 

-P 12 = 

Pll ~ P22 + Pl2 2 » 

Pl2( 1) = 

-P 21 = 

Pll + Pll P22 + 1 

, p 21 ( i) = i 

-P 22 = 

P21 + Pl2 + Pl2 P22 

> P 22 ( i) = 


Equations (144) and (145) correspond to a single fourth-order equation in 
xj for which a solution can also be obtained by the classical Laplace trans- 
form method. The procedure is straightforward; however, evaluating the 
residues at the poles s,., i=l, 2, 3, 4 is very tedious since the initial values 
of x 1 and jq are not specified. A computer algorithm for evaluating the 
inverse transform would be helpful. Unfortunately, the solution does not 
exist on the infinite interval [0, °°] since the system possesses poles in the 
right-half complex plan 
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CHAPTERS 


CONCLUSIONS AND RECOMMENDATIONS 


Few classification schemes are perfect since there is usually overlap, 
and membership in a single category is often fuzzy. However, the scheme 
given in this report appears to be logical and easy to apply. Most of the well- 
known computational methods for optimal control belong to one of the categories, 
including steepest descent, quasilinearization, boundary condition iteration, 
and the "direct" second variation method. The classification provides geometric 
insight for the design of new algorithms, as evidenced by the second variation 
method for Approximation in X . 

The conjugate direction algorithm for Approximation in X" requires a 
moderate amount of instructions and computation. The method inherits all 
the advantages of the particular quasi-Newton method which is employed, 
including rapid and stable descent to a minimum. In contrast, methods 
involving the Newton-Raphson method may suffer from ill conditioning of the 
Jacobian matrix and unstable descent, in addition to the difficulty in estimating 
second partial derivatives of the cost functional by divided differences. 

The monotone convergence algorithm for Approximation in U requires 
a modest amount of instructions and computation, even for high-order systems. 
Numerical results indicate that the algorithm is preferable to existing techniques, 
including Runge-Kutta integration of the Riccati equations or the "Automatic 
Synthesis Program" (ASP). In the latter two methods, convergence is linear, 
whereas in the new method convergence is both monotonic and quadratic. 

The second variation algorithm for Approximation in X requires a 
large amount of instructions and computation. However, the transition matrix 
approach offers the advantage that only functions evaluated at the terminal time 
are saved between control updates. For problems in which the Euclidean 
dimension of the control space U is less than that of the state space X , 
storage requirements are less than former second variation methods. Solution 
of the Accessory Problem by the sweep method, on the other hand, requires 
nearly the same amount of computation as former techniques. 

Future research in the area of this report should concern extending 
the algorithms for Approximation in X and Approximation in U to handle 
terminal equality constraints on the state with unspecified terminal time. For 
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the second variation method, it is expected that a slight redefinition of the 
Hilbert space inner product as in Reference 10 would be needed for the 
latter purpose. Further consequences of the inverse mapping should be 
investigated for high-order time-varying systems by application of the 
algorithm to physical problems. 

Further research should involve the application of the algorithms to 
constrained nonlinear problems. At present there are a number of methods 
for handling constraints that restrict elements of the spaces U and X to 
closed subsets of Hilbert space [44], Among these are (1) methods of 
feasible directions such as the "projected gradient scheme" ; (2) penalty 
function methods; (3) methods of set approximation such as the Ritz method: 
(4) duality methods such as the "moment methods" based on the Krein 
L-problem and also methods based on the Kuhn Tucker conditions; (5) meth- 
ods of optimal evolution including the well-known Dynamic Programming 
technique. Methods (3) and (4) appear to be most promising for future 
development. 
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APPENDIX A 


SOLUTION OF THE LINEAR EQUATION 

_p(k) . p <k) A (k) + ( A (k)J p (k) + q(IO 

In this section we consider the computational solution of the linear 
differential equation that was obtained by application of the method for 
monotone Approximation in U . 

Suppose A, B, Q, R , are constant matrices, and t t — + °° . 

(k) 

Consider a translation of the time origin such that t 0 = 0 . Since A' is a 
stability matrix (that is, the system equation (67) with the control function 

• (k) 

defined by equation (71) is asymptotically stable), lim P' ' (t) = © , and 

t — - 

we may obtain the optimal control approximation from equation (71) by 
solving the linear matrix equation 

P (k) A (k) + (a ( 10 )* P (k) +Q (k) = 0 , (A. 1) 

(k) (k) 

where A x ' and Q' ' are given by equations (73) and (74), and 


P <k> = lim p <k) (t) 


Then equation (A. 1) may be written as fn(nfl) linear algebraic equations 
whose solution is 


where 


(k) 


(k) 


■< 


-1 


D 


-(■ 


,< k n 

q <« 



(k) 

D (k) - 

o (k) 

(k) 

ii *• 

‘ ‘ ’ P in ’ 

P 22 ’ * * 

' ’ P 2n ” 

(k) 

„<k) 

o (k) 

a (k) . 

ii ’* 

* • ’ q in ’ 

q 22 ’• • 

‘ ’ q 2n ’* 



(A. 2) 
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and is an |-n(n+l) x£n(n+l) constant matrix whose elements are 


linear combinations of elements of A 
is as follows: 


<k) 


For example, if n = 3 , D 


00 


,(k) 


, ( k> 

2a ll 

2a (k) 

2a 21 

2a (k) 
2 31 


0 

0 

( k I 
a i2 

(k) (k) 

a t- a 

11 22 

a (k) 

21 


(k) 

a 21 

(k) 

a 31 

(k> 

“13 

a (k) 

23 

a (k) * 
11 

a <k) 

33 

0 


0 

o < k) 
2a !2 

0 


2a (k) 

2a 22 

“S’ 

0 

, <k) 

d 13 

<k) 

a i2 


(k) 

a 23 

(k) 

a 22 

0 

0 

2a (k) 
2 13 


0 

2a (k) 

23 


(k) 

*33 


(k) 

31 



For the finite interval case, t x - t 0 < «= and the elements of A 

,00 


(k) 


and Q are functions on T . Let us consider the following expansions: 


P (k) (t) = Yj P^ k) (tj - t) n 
n-1 

A (k) (t) = ZA< k) (t,-t) n 
n=0 

Q (k) (t) = Y Q <k) (tj - t) n 
n=0 n 


(A. 3) 


(A. 4) 


(A. 5) 


( k ) ( lc ) ( k ) 

where P .A , and Q are constant n \ n matrices. Substituting 

" n n (k) (k) (k) 

the latter expressions for P v , A v , and Q into equation (72) 

and equating coefficients of (t| - t)” . n - 0. 1. 2, .... we obtain the 
following set of equations: 
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Equation (A. 6) may be solved recursively for each value of k and provides 
one means for solving the linear matrix differential equation (72). Motivated 
by this result, a similar approach was used to generate simple optimal test 
problems and is described in Appendix C. 

In the general case of time-varying A, B, Q, or R , equation (72) 
must be integrated numerically using a digital computer. Note that equation 

(k) 

(72) is integrated backward from the final condition P' (t t ) = 0 . 
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APPENDIX B 


SOLUTION OF THE ACCESSORY PROBLEM 


For the computational solution of the boundary value problem, equation 
(103) through (111) , it is convenient to consider the following two cases: 

(1) h > 0 ; (2) h ^ 0 and L > 0 . When applicable, the procedure 

described for Case 1 is easier to implement and requires less computation. 


Case 1. h > 0 „ 

TRANSITION MATRIX APPROACH 

2n 

Let us rewrite equation (103) as the following IR -valued system of 
equations: 


5i - S6z + w 

where 


Sz 




(B. 1) 


and subject to the boundary conditions, equations (108) and (109), 

5x(t 0 ) = 6 2) 

6x(tj) =h _1 (t 1 ) N(tj) 6x(tj) + h -1 (t, ) r(tj) . (B.3) 

The solution to equation (B. 1) can be written in terms of the initial conditions 
6x(t 0 ) and <5x(t 0 ) as 

<5x(t) =4> 11 (t,t 0 ) <5x(t 0 ) +$ 12 (t,t 0 ) 6x(t 0 ) +6x(t) , (B.4) 

<5x(t) =$ 21 (t,t 0 ) 5x(t 0 ) +<f> 22 (t,t 0 ) Sx(t 0 ) +6x(t) , . (B. 5) 
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where <i> 1J (t, are n x n submatrices of the transition matrix for the linear 
system [eq. (B. 1)], 6x(t 0 ) is unknown, and Sx , 6x , represent the unknown 
"forced solutions" of equation (B. 1). Evaluating equation (B. 5) at t = tj and 
using equations (B. 2) and (B.3), we obtain 

$ 22 (ti » to) 6x(t 0 ) = 6x(tj) - dxltj) 

= h -1 !^) N(tj) fodtj) + h _1 (tj ) r(tj) - dxltj) 

(B. 6) 

= h _1 (tj) Nltj) # 12 (ti , t 0 ) <5x(t 0 ) 

+ h~ I (t 1 ) N(tj ) dxltj ) - 6x(tj ) +h -1 (tj) rlt^) 

Equation (B. 6) may be solved for the "missing" initial condition, 

6x(t 0 ) = ( h(tj) $ 22 (tt , to) -N(tj) * 12 <t, , tg))" 1 x 

(N(tj) 6x(tj) - 6x(t[) +r(t t )) . (B.7) 


If we can determine 6x(tj) , 6x(tj) and $^(ti , t 0 ) , i = 1,2 , the 
missing initial condition &x(t 0 ) is specified by equation (B. 7). A procedure 
for computing these unknown functions is as follows. Let fix ( t 0 ) = fix ( t 0 ) = 6 
inequation (B. 1). Then 6x(t) = 6x(t) and <5x(t) =6x(t) . By integrating 
equation (B. 1) with these initial conditions until t = t* , we obtain the "forced 

solution" 6x(tj) . To compute the submatrices <& (t t , t 0 ) , i = 1, 2, let 

k = 9 in equation (B. i) so that the "forced solutions in equations (B. 4) and 

(B. 5) are zero. Then by integrating equation (B. 1) with the initial conditions 

6x(t 0 ) = 6 and x(t 0 ) = {A . , A , . . . , A } , where A is defined by 

lj 2j nj lj 


A.. = 
h) 


i = 3 
i * 3 


the jth column of $(tj , t 0 ) is obtained. By repeating this procedure for 
j = i,2,...,n,we obtain the last n columns of $(ti , t 0 ) as required. 

A total of 2n(m + 2) differential equations must be integrated where 
m = rank h . Existence of the inverse in equation (B. 7) is insured by the 
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absence of conjugate points in the interval T [45.46] and, together with the 
strong Legendre condition h > 0 . constitutes sufficient conditions for a solu- 
tion to the Accessory Problem. 

THE SWEEP METHOD 

Alternatively, boundary value problems of this type may be solved 
by the sweep method [6j. Consider the inhomogeneous Riccati transformation 


<5x = P, 5x + qj 


(B. 8) 


where Pj is an n x n matrix, and q, is an IR n -valued function. It is 
sufficient to assume that elements of P[ and qj are real-valued continuous 
functions on T . Then 

<5x = P t 5x + Pj 6x + qj 


= (Pj + F?) <5x + (P^ + qj) 


= h _1 (F<5x + G5x + k) 

In the last step we have used equation (103). Thus, the following equations for 
the "backward sweep" are obtained: 

-hPj - hPj 2 - FPj - G , h(t,) Pi(tj) =N(t,) (B. 9) 

-hq,= (hP t - F) q,- k , h ( t t ) q 4 ( tj ) = r ( tj ) . (B. 10.) 

During the "forward sweep" we must integrate 

ox Pjdx+qj , 6x(t 0 ) -6 (B.ll) 

to obtain ox cm T . From equations (94) and (B. 8) the control increment 
may be expressed in terms of the state increment as follows: 

ou 5' ox + z, 


(B. 12) 


where the "feedback" and "feedforward" gains are. respectively, 


Y,=g (Pj'f) ( B - 13 ) 

fc A 

z, = g e qj (B. 14) 

Existence of a solution to equation (117) is insured by the conjugate 
point condition, [6,45]. In contrast to Riccati equations encountered in the sec- 
tions, State Regulator Problems, Chapter 2, and Methods Based on the Second 
Variation, Chapter 3, Pj in equation (B. 9) is not necessarily symmetric, and 
hence n(n+2) equations must be integrated to obtain 6u . In comparison, 
the Riccati transformation for the second variation method described in the 
section. Descent Algorithms for Hilbert Space, Chapter 3, results infn(n+9) 
equations 


Case 2. h > 0 and L uu > 0. 


In this case equations (103) through (11) may be written as the 
valued system of linear equations 


IR 


2n 



where 6x = (6x^ , . . . , 6x_) , 6x - (^ x n+ ^ , . . . . 6x ) satisfy 

6x = 6x , i = 1. 2. . . . 2n - 1 with boundary conditions 

i + 1 i ' 


<5x a (t 0 ) = Q 


(B. 16) 


C(t t ) 6x a (t,) = D(t 1 ) 6x b (t 1 ) - e ( t| ) (B. 17) 

Notice that by definition <5x a ■= 6x and 6x* J -- 6x . Matrices , C , D 
a b 

and d . d are n x n and n x 1 , respectively, with elements that are 
functions on T . It is assumed that D(t t ) is invertible. 
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TRANSITION MATRIX APPROACH 
In analogy to Case 1 let 


/ 6x a (t 1 ) 
\ 6x b (t 1 ) 



^Mtj.to) $ 12 (t! , to)\ / 6x (to) 


# 21 (ti,t 0 ) # 22 (tj , t 0 )/ \5x (t 0 ) 


(t^ , s) $ 12 (tj , s) 


t n V 21 ^ > S ) ^ 22(t l » t ( 



(B. 18) 


where ^(tj , t 0 ) are n x n submatrices of the transition matrix for the 
linear system equation (B.15). Combining equations (B. 18), (B. 16), and 
(B. 17) , we obtain 


C(tj) ( * 12 (tt , to) 6x b (t 0 ) + f , s) d a (s) + -h 12 ^ , s) d b (s)) ds 

V to 


= D^)^# 22 ^ , t 0 ) 5x b (t 0 ) + f ^ 21 (t! , s) d a (s) + $ 22 (t, , s) d b (s)) dsj - 


e(t) . 
(B. 19) 


Equation (B. 19) may be solved for the missing initial condition as 
x b (t 0 ) = ( D(tj) dh 22 (tq , t 0 ) - C (tj ) * 12 (t! , to ))" 1 < 

C (tj ) S 11 (tj , s) - D(tj) d’ 21 !^ , s))d a (s) 

+(c(ti) ^^(tj , s) - D(tj) d> 22 (t[ , s))d b (s) 

(B. 20) 
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The linear system equation (B. 15) , together with initial conditions in 
equations (B. 16) and (B. 17), constitutes an initial value problem, which, 

ci b 

if the inverse exists, may be integrated forward to obtain 6x and <5x on T . 

Existence of the inverse in equation (B. 20) is insured by the absence of 
conjugate points in the interval T [45,46] and, together with the strong 
Legendre condition L > 0 , represents sufficient conditions for a solution 

to the Accessory Problem. 

SWEEP METHOD 

Consider the following inhomogeneous Riccati transformation: 

<5x b = P 2 5x a + q 2 , (B. 21) 

n 

where P 2 is an n x n matrix on T , and q 2 is an IR -valued function 
on T . Then 

_ .b -jL , a a . 

<5x = P 2 <5x + P 2 <5 x + q 2 

= P 2 6x a + P 2 A n 5x a + A i2 P 2 Sx a + q 2 + d & + q 2 (B. 22) 

= A 21 Sx a +A 22 P 2 <5x a + q 2 + d & , 

where we have used equations (B. 21) and (B. 15). Thus the following equations 
for the "backward sweep" are obtained: 

P 2 = - P 2 A 11 +A 22 P 2 - P 2 A 12 P 2 + A 21 , DU}) P 2 (t 1 ) = C (tj ) (B. 23) 

q 2 = (A 22 - P 2 A 12 )q 2 - P 2 d a + d b , D(tj) q 2 (t 1 ) = eltj) . (B. 24) 

During the ’’forward sweep” we must solve 

6x a = (A 11 +A 12 P 2 ) 6x a +A 12 q 2 +d a , Sx a (t 0 ) = 6 . (B.25) 
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H 3, 

Having obtained the IR -valued function 6x , the control increment is 

5u = Y 2 5x + z 2 , (B. 26) 

where the "feedback" and "feedforward" gains are, respectively, 

Y 2 = g (A 11 +A 12 P 2 - f ) (B. 27) 

t X. 

z 2 = g £ A 12 q 2 . (B. 28) 

Existence of a solution to equation (B. 23) is insured by the conjugate point 
condition [45, 6] . Since A 12 and A 21 are not necessarily symmetric, 

Si 

n(n+2) equations must be integrated to obtain &x = <5x , 
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APPENDIX C 


A CLASS OF STATE REGULATOR PROBLEMS 
FOR TESTING COMPUTATIONAL ALGORITHMS 


This section describes a class of simple first-order time-varying 
dynamical systems and the optimal solution with respect to quadratic cost as 
specified by the Riccati solution of the section, State Regulator Problem, 
Chapter 2. These problems have been designed to provide a class of simple 
problems with nontrivial closed-form solutions for the testing of computational 
algorithms and are based on Reference 47. 

The system dynamics are assumed to be modeled by a linear first- 
order differential equation, 


x = Ax + u , (C. 1 ) 

where A is a continuous function on [0 , tj] that will be determined; 

x eLf (T) and ueL^fT) are the state and control functions, respectively. 
The cost functional is 


J(u , x) = i J (Qx 2 +u 2 ) dt , (C. 2) 

0 

where Q is a continuous function 1 0 , tj which will also be determined. 
From the section, State Regulator Problem, the optimal feedback solution to 
equations (A.l) and (A. 2) is 

u = - Px , (C. 3) 

where P satisfies the Riccati differential equation 

-P = 2PA - P 2 +Q (C. 4) 

with boundary condition P(tj) = 0 . Equation (C.4) must be solved in 
reverse time starting with the boundary condition and integrating backwards. 
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(C. 5) 


Example 1 

Let a , q p be constants and suppose that 
n n n 

Ait) = a 0 + aj (tj - t) 

2 

Q(t) = E q_(ti - t) n (c.6> 

n=0 

P(t) = pj (tj - t) (C. 7) 

By substituting equations (C. 5) and (C.7) into equation (C.4) and collecting 

coefficients of (tj - t) n , n = 0, 1, . . . we obtain the following constraints 
on our choice of system parameters: 


Pi = To 


(C. 8) 

Ti + 2a 0 Pi 

= 0 

(C. 9) 

q 2 + 2ai Pi 

- P? = o 

(C. 10) 


Hence, given the system equations (C. 1) and (C.5), equations (C.2) and 
(C. 6) specify the cost functional for linear feedback. For example, suppose 
that 


A(t) - (-1 +$(t, - t)) (C. 11) 

Then a 0 - 1 , \ . If pj = \ 5 Qo = 2 » then from equations 

(C. 9) and (C. 10) we find that the controlled system is specified by 


X = 

(- 1 -J(tj - t) ) X 

(C. 12) 

.A 

J (u . 

tl 

x) = \ f { \ + (tj - t) + l (t, - t) 2 ) x 2 dt ^x 2 (0) 

(C. 13) 


0 


A 

U - 

- |(tj - t)x 

(C.14) 



Example 2 


Let a^ , q^ , be constants and suppose that 

, a(ti - t) 

A(t) = a 0 + aje n 

cult - t) 2a (L - t) 

Q(t) = q 0 + q t e 1 + q 2 e 1 


P(t) = p 0 + p t e 


a 


(ti - 1> 


(C* 15) 
(C. 16) 

(C. 17) 


Substituting equations (C. 15) through (C.17) into equation (C.4) and collecting 

coefficients of e na ^ ^ , n = 0 , 1 , . . . we obtain the following 

constraints 


q 2 + 2p!aj - p? = 0 

(C. 18) 

apt = q t + 2a 0 Pi + 2a! p 0 - 2p 0 P! 

(C. 19) 

q 0 + 2a 0 p 0 - p§ = 0 

(C. 20) 


In addition, since P(t t ) = 0 , 


Po + Pi - 0 . (C. 21) 

Equations (C. 18) through (C.21) specify the requirements on our choice of 
cost functional equation (C. 16) for the system equation (C. 15) and control 
in terms of equation (C. 17). For a < 0 and tj — + °° the optimal control 
will approach that for the constant coefficient system 


x = a 0 x + u 


(C . 22) 


J(u,x) = f (q 0 x 2 +u 2 )dt 
0 


(C. 23) 
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Example 3 

Let a^, q^, be constants and suppose that 

A Oft 

A = a 0 + a^e 
_ cet 2o!t 

Q = q 0 + Qi e +q2 e 

_ at 

P = Po + Pi e 

Substituting in equation (C.4) and collecting coefficients as above, we 

2a 0 Po + q 0 - Po = 0 

- a pj = 2a! Po + 2aoPj + q t - 2a 0 pj 
2^ + q 2 - Pi = 0 


OlU 

Po + Pie n = o 

Consider the case of infinite time interval [0 , °°] and a = -1 
Moreover, suppose that 

a 0 = 3 q 0 = o qi = i Po = o 

a i - s q2 = o Pi =1 

Then the optimal solution is characterized by 


x = + e fc ) x + « 

Q ( t) — 

P(t) = e _t 

a -t 

u = - e x 


{C. 24) 
(C. 25) 

(C. 26) 
obtain 

(C. 27) 

(C. 28) 

(C. 29) 

(C. 30) 


(C. 31) 
(C. 32) 
(C. 33) 
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(C. 34) 



I 


It is interesting that the dynamics of the closed-loop system are 
x = - e fc )x 

an optimal system which becomes unstable within one second! Another 
example for this case is generated by the following choice: 

a o = i q o = 0 q i = | Po = 0 

% = i q 2 = 0 Pi = 1 

for which the optimal system is characterized by 

x = i (~i + e -t )x + u 
Q(t) = | e t 

P(t) = e _t 

A -t 

u = - e x 

Example 4 

Let 

A ( t) = a 0 + aj sin ait + a 2 cos cvt 

Q(t) = q 0 + q i sin cot + q 2 cos cot + qg sin 2cot + cos 2cot 

P(t) = p 0 + Pi sin cot + pg cos cot 

The resulting constraints are obtained as follows: 

2aflP 0 + q 0 - p 5 a,pj + a 2 P 2 - i(p? - p!) = 0 


I 


(C. 35) . 


(C. 36) 
(C. 37) 

(C. 38) 
(C. 39) 

(C. 40) 
(C. 4i) 
(C. 42) 

(C. 43) 
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wp 2 + 2(a 0 pj + aiPo - p 0 Pi) + q t = 0 

( C . 44) 

-«Pi + 2(aoP2 + a 2 P 0 - p 0 p 2 ) + q 2 = 0 

(C. 45) 

aiP 2 + 32 Pi + <k + P 1 P 2 = 0 

(C. 46) 

-a 1 pi + a2P2 + q4 + £(pj - p|) = 0 

(C. 47) 

Po + pj sin wtj + P 2 cos cutj = 0 

(C.48) 


Equation (C.43) through (C.48) may be used to generate optimal test 
problems, as in the previous example. 

Example 5 

Let a^, q^, p^ be constants and suppose that 

A ( t) = fj a e 11 " 4 (C. 49) 

n=0 n 


_ /A . v n ol t 

Q(t)= Ij a e (C. 50) 

n=0 


V nat 

P(t)= h P e (C. 51) 

n=0 


The constraints on a^ and q^ are obtained by substituting equations (C. 49) 
through (C.51) into equation (C.4). Thus 


Po = 0 


■ ap 


n+1 


1 

n+1 



n 


z 

v=i 


(2p a 
' v n-v 


P P 

v n-v , 


n - 0, 1, 2, 


(C. 52) 
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where the terminal constraint p 0 = 0 has been combined with equation (C. 51) 
to yield 


P(t) = 


OO / 

2 ( p „( eI1 

n=l \ 


a t 


- e 


no; tj 


Example 6 

Let a , q , p be constants and suppose that 
n n n 


A(t) = J at” 

« n 
n=0 

OO 

Q<t>= Z v n 

n=0 

OO 

p(t)= Z p n (ti - t) n 

n=l 


In a similar manner we obtain constraints on a and q in the form 
recursion relations: 

Po = 0 


-ap 


n+1 



n 


z 


1^=1 


(2p a 
\ v n-v 


~ V 


n-u 



n = 0, 1, 


which may be used to generate optimal solutions to the problem. 


George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama, 35812, Sept. 15, 1970 

125-17-05-00-62 


(C. 53) 

(C. 54) 

(C. 55) 

(C. 56) 
of the 

2 , • . . 

(C. 57) 
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